P2: Análisis de Supervivencia
Naira Maria Chiclana García
Noviembre 2018
1.Objetivos
Obtendremos los datos reales de pacientes con cancer de mama de The Cancer Genome Atlas y prepararemos estos(convertirlos e imputarlos) para su consecuente análisis. Una vez los datos tengan el formato correcto, los estratificaremos de la forma en la que suelen tratarse en el ámbito médico. Por grupos tumorales, Luminal, grupos nodales, estado nodal y grupos de edad. 2
Haremos un barrido rápido por algunas las opciones de análisis (estadístico 3.1, paramétrico 3.3 y no paramétrico 3.2) que podemos aplicar en análisis de supervivencia. Para ello usaremos principalmente el paquete
survminer.Cuando tengamos el dataset final adecuado para trabajar con los datos correctos y las estratificaciones incluidas pasaremos a ver que variables nos interesan y cuales no, es decir, cuales son útiles para explicar los riesgos y la supervivencia. Este proceso lo llevaremos a cabo gracias al análisis multivariante del modelo semiparamétrico Cox. Realizaremos una búsqueda exhaustiva de las mejores covariantes (de forma univariante y multivariante), y después de descartar aquellos modelos que no cumplan la proporcionalidad de riesgos, nos quedaremos con aquella combinación de variables que explique mejor el riesgo. 3.4
Una vez tengamos las variables con las que trabajaremos más en detalle, haremos un análisis exploratorio visual de como interactuan y explican el riesgo diferentes combinaciones de variables 3.5. De aquí sacaremos conclusiones de las cuales le haremos la consiguiente prueba de hipótesis 3.6.
Analizaremos más en detalles los supuestos aceptados con más significancia. 3.7
2.Preparación datos
Instalación de paquetes y carga de librerias:
install.packages("survival")
install.packages("KMsurv")
install.packages("survMisc")
install.packages("survminer")
install.packages("UCSCXenaTools")
install.packages("TCGAretriever")library(KMsurv)
library(survMisc)
library(survminer)
library(flexsurv)
library(dplyr)
library(UCSCXenaTools)
library(TCGAretriever)
library(ggfortify)Los extraeremos de TCGA, y después de imputarlos, los pondremos en el formato más conveniente para trabajarlos.
2.1. Importar los datos de TCGA
Obtendremos los datos de TCGA The Cancer Genome Atlas, usando el paquete TCGAretriever.
all_studies<-get_cancer_studies() #identifier, tittle and description of the study
str(all_studies)## 'data.frame': 233 obs. of 3 variables:
## $ cancer_study_id: chr "paac_jhu_2014" "all_stjude_2016" "laml_tcga_pub" "laml_tcga_pan_can_atlas_2018" ...
## $ name : chr "Acinar Cell Carcinoma of the Pancreas (Johns Hopkins, J Pathol 2014)" "Acute Lymphoblastic Leukemia (St Jude, Nat Genet 2016)" "Acute Myeloid Leukemia (TCGA, NEJM 2013)" "Acute Myeloid Leukemia (TCGA, PanCancer Atlas)" ...
## $ description : chr "Mutation data from whole exome sequencing of 23 surgically resected pancreatic carcinomas with acinar differentiation." "Whole-genome and/or whole-exome sequencing was performed for ERG-altered B-ALL cases." "TCGA Acute Myeloid Leukemia, analysis of 200 adult cases; raw data at the <A HREF=\"https://tcga-data.nci.nih.g"| __truncated__ "Comprehensive TCGA PanCanAtlas (https://gdc.cancer.gov/about-data/publications/pancanatlas) data from 11k cases"| __truncated__ ...head(all_studies[,2],20) #tittles of first 20 studies## [1] "Acinar Cell Carcinoma of the Pancreas (Johns Hopkins, J Pathol 2014)"
## [2] "Acute Lymphoblastic Leukemia (St Jude, Nat Genet 2016)"
## [3] "Acute Myeloid Leukemia (TCGA, NEJM 2013)"
## [4] "Acute Myeloid Leukemia (TCGA, PanCancer Atlas)"
## [5] "Acute Myeloid Leukemia (TCGA, Provisional)"
## [6] "Adenoid Cystic Carcinoma (FMI, Am J Surg Pathl. 2014)"
## [7] "Adenoid Cystic Carcinoma (MDA, Clin Cancer Res 2015)"
## [8] "Adenoid Cystic Carcinoma (MSKCC, Nat Genet 2013)"
## [9] "Adenoid Cystic Carcinoma (Sanger/MDA, JCI 2013)"
## [10] "Adenoid Cystic Carcinoma of the Breast (MSKCC, J Pathol. 2015)"
## [11] "Adrenocortical Carcinoma (TCGA, PanCancer Atlas)"
## [12] "Adrenocortical Carcinoma (TCGA, Provisional)"
## [13] "Ampullary Carcinoma (Baylor College of Medicine, Cell Reports 2016)"
## [14] "Bladder Cancer (MSKCC, Eur Urol 2014)"
## [15] "Bladder Cancer (MSKCC, JCO 2013)"
## [16] "Bladder Cancer (TCGA, Cell 2017)"
## [17] "Bladder Cancer, Plasmacytoid Variant (MSKCC, Nat Genet 2016)"
## [18] "Bladder Urothelial Carcinoma (BGI, Nat Genet 2013)"
## [19] "Bladder Urothelial Carcinoma (Dana Farber & MSKCC, Cancer Discov 2014)"
## [20] "Bladder Urothelial Carcinoma (TCGA, Nature 2014)"head(all_studies[,1],20) #id of first 20 studies## [1] "paac_jhu_2014" "all_stjude_2016"
## [3] "laml_tcga_pub" "laml_tcga_pan_can_atlas_2018"
## [5] "laml_tcga" "acyc_fmi_2014"
## [7] "acyc_mda_2015" "acyc_mskcc_2013"
## [9] "acyc_sanger_2013" "acbc_mskcc_2015"
## [11] "acc_tcga_pan_can_atlas_2018" "acc_tcga"
## [13] "ampca_bcm_2016" "blca_mskcc_solit_2014"
## [15] "blca_mskcc_solit_2012" "blca_tcga_pub_2017"
## [17] "blca_plasmacytoid_mskcc_2016" "blca_bgi"
## [19] "blca_dfarber_mskcc_2014" "blca_tcga_pub"#filter rows from all_studies which name contains breast
df.breast<-dplyr::select(filter(all_studies, grepl('Breast', all_studies$name)),everything())
df.breast$cancer_study_id## [1] "acbc_mskcc_2015" "brca_metabric"
## [3] "breast_msk_2018" "bfn_duke_nus_2015"
## [5] "brca_bccrc" "brca_broad"
## [7] "brca_sanger" "brca_tcga_pub2015"
## [9] "brca_tcga_pub" "brca_tcga_pan_can_atlas_2018"
## [11] "brca_tcga" "brca_bccrc_xenograft_2014"
## [13] "brca_mbcproject_wagle_2017"#from all the studies, we choose brca_tcga->brca_tcga_all
all_case_lists<-get_case_lists("brca_tcga")
colnames(all_case_lists)## [1] "case_list_id" "case_list_name" "case_list_description"
## [4] "cancer_study_id" "case_ids"all_case_lists[,1:2] #id and name## case_list_id
## 1 brca_tcga_3way_complete
## 2 brca_tcga_sequenced
## 3 brca_tcga_methylation_all
## 4 brca_tcga_all
## 5 brca_tcga_protein_quantification
## 6 brca_tcga_cna
## 7 brca_tcga_methylation_hm27
## 8 brca_tcga_methylation_hm450
## 9 brca_tcga_mrna
## 10 brca_tcga_rna_seq_v2_mrna
## 11 brca_tcga_rppa
## 12 brca_tcga_cnaseq
## case_list_name
## 1 All Complete Tumors
## 2 All Sequenced Tumors
## 3 All tumor samples with methylation data
## 4 All Tumors
## 5 Protein Quantification (Mass Spec)
## 6 Tumor Samples with CNA data
## 7 Tumor Samples with methylation data (HM27)
## 8 Tumor Samples with methylation data (HM450)
## 9 Tumor Samples with mRNA data (Agilent microarray)
## 10 Tumor Samples with mRNA data (RNA Seq V2)
## 11 Tumor Samples with RPPA data
## 12 Tumor Samples with sequencing and CNA datadf.tcga<-get_clinical_data("brca_tcga_all")
str(df.tcga)## 'data.frame': 1105 obs. of 110 variables:
## $ CASE_ID : chr "TCGA-A7-A3J0-01" "TCGA-OL-A66N-01" "TCGA-AQ-A0Y5-01" "TCGA-E9-A22H-01" ...
## $ AGE : chr "62" "59" "70" "42" ...
## $ AJCC_METASTASIS_PATHOLOGIC_PM : chr "M0" "MX" "MX" "M0" ...
## $ AJCC_NODES_PATHOLOGIC_PN : chr "N0" "N3" "N2a" "N1" ...
## $ AJCC_PATHOLOGIC_TUMOR_STAGE : chr "Stage IIA" "Stage IIIC" "Stage IIIA" "Stage IIB" ...
## $ AJCC_STAGING_EDITION : chr "7th" "7th" "7th" "7th" ...
## $ AJCC_TUMOR_PATHOLOGIC_PT : chr "T2" "T3" "T2" "T2" ...
## $ BRACHYTHERAPY_TOTAL_DOSE_POINT_A : chr "" "" "" "" ...
## $ CANCER_TYPE : chr "Breast Cancer" "Breast Cancer" "Breast Cancer" "Breast Cancer" ...
## $ CANCER_TYPE_DETAILED : chr "Breast Invasive Mixed Mucinous Carcinoma" "Breast Invasive Lobular Carcinoma" "Breast Invasive Ductal Carcinoma" "Breast Invasive Ductal Carcinoma" ...
## $ CENT17_COPY_NUMBER : chr "" "" "" "" ...
## $ DAYS_TO_BIRTH : chr "-22672" "-21608" "-25793" "-15504" ...
## $ DAYS_TO_COLLECTION : chr "78" "505" "212" "38" ...
## $ DAYS_TO_DEATH : chr "" "" "172" "" ...
## $ DAYS_TO_INITIAL_PATHOLOGIC_DIAGNOSIS : chr "0" "0" "0" "0" ...
## $ DAYS_TO_LAST_FOLLOWUP : chr "76" "413" "" "0" ...
## $ DFS_MONTHS : chr "10.28" "26.02" "" "40.47" ...
## $ DFS_STATUS : chr "DiseaseFree" "DiseaseFree" "" "DiseaseFree" ...
## $ DISEASE_CODE : chr "" "" "" "" ...
## $ ER_POSITIVITY_SCALE_OTHER : chr "" "" "" "" ...
## $ ER_POSITIVITY_SCALE_USED : chr "" "" "" "" ...
## $ ER_STATUS_BY_IHC : chr "Positive" "Positive" "Positive" "Positive" ...
## $ ER_STATUS_IHC_PERCENT_POSITIVE : chr "90-99%" "80-89%" "70-79%" "<10%" ...
## $ ETHNICITY : chr "NOT HISPANIC OR LATINO" "NOT HISPANIC OR LATINO" "NOT HISPANIC OR LATINO" "NOT HISPANIC OR LATINO" ...
## $ FIRST_SURGICAL_PROCEDURE_OTHER : chr "" "" "" "" ...
## $ FORM_COMPLETION_DATE : chr "3/13/12" "6/11/13" "4/7/11" "8/2/11" ...
## $ FRACTION_GENOME_ALTERED : chr "0.224023880" "0.209382668" "0.496072001" "0.103942673" ...
## $ HER2_AND_CENT17_CELLS_COUNT : chr "" "" "" "" ...
## $ HER2_AND_CENT17_SCALE_OTHER : chr "" "" "" "" ...
## $ HER2_CENT17_RATIO : chr "" "1.1" "3.8" "" ...
## $ HER2_COPY_NUMBER : chr "" "" "" "" ...
## $ HER2_FISH_METHOD : chr "" "" "" "" ...
## $ HER2_FISH_STATUS : chr "" "Negative" "Positive" "" ...
## $ HER2_IHC_PERCENT_POSITIVE : chr "" "" "" "<10%" ...
## $ HER2_IHC_SCORE : chr "1" "" "2" "" ...
## $ HER2_POSITIVITY_METHOD_TEXT : chr "" "" "" "" ...
## $ HER2_POSITIVITY_SCALE_OTHER : chr "" "" "" "" ...
## $ HISTOLOGICAL_DIAGNOSIS : chr "Mucinous Carcinoma" "Infiltrating Lobular Carcinoma" "Infiltrating Ductal Carcinoma" "Infiltrating Ductal Carcinoma" ...
## $ HISTOLOGICAL_SUBTYPE : chr "" "" "" "" ...
## $ HISTORY_NEOADJUVANT_TRTYN : chr "No" "No" "Yes" "No" ...
## $ HISTORY_OTHER_MALIGNANCY : chr "No" "No" "Yes" "No" ...
## $ ICD_10 : chr "C50.9" "C50.9" "C50.9" "C50.9" ...
## $ ICD_O_3_HISTOLOGY : chr "8480/3" "8520/3" "8500/3" "8500/3" ...
## $ ICD_O_3_SITE : chr "C50.9" "C50.9" "C50.9" "C50.9" ...
## $ IHC_HER2 : chr "Negative" "" "Positive" "Positive" ...
## $ IHC_SCORE : chr "" "" "" "" ...
## $ INFORMED_CONSENT_VERIFIED : chr "YES" "YES" "YES" "YES" ...
## $ INITIAL_PATHOLOGIC_DX_YEAR : chr "2011" "2011" "2010" "2011" ...
## $ INITIAL_WEIGHT : chr "110" "310" "160" "500" ...
## $ IS_FFPE : chr "NO" "NO" "NO" "NO" ...
## $ LYMPH_NODES_EXAMINED : chr "YES" "YES" "YES" "YES" ...
## $ LYMPH_NODES_EXAMINED_HE_COUNT : chr "0" "13" "5" "1" ...
## $ LYMPH_NODES_EXAMINED_IHC_COUNT : chr "" "" "" "0" ...
## $ LYMPH_NODE_EXAMINED_COUNT : chr "2" "15" "17" "8" ...
## $ MARGIN_STATUS_REEXCISION : chr "" "" "" "" ...
## $ MENOPAUSE_STATUS : chr "Post (prior bilateral ovariectomy OR >12 mo since LMP with no prior hysterectomy)" "Post (prior bilateral ovariectomy OR >12 mo since LMP with no prior hysterectomy)" "Post (prior bilateral ovariectomy OR >12 mo since LMP with no prior hysterectomy)" "Pre (<6 months since LMP AND no prior bilateral ovariectomy AND not on estrogen replacement)" ...
## $ METASTATIC_SITE : chr "" "" "" "" ...
## $ METASTATIC_SITE_OTHER : chr "" "" "" "" ...
## $ METASTATIC_TUMOR_INDICATOR : chr "NO" "" "NO" "NO" ...
## $ METHOD_OF_INITIAL_SAMPLE_PROCUREMENT : chr "Core needle biopsy" "Core needle biopsy" "Core needle biopsy" "Tumor resection" ...
## $ METHOD_OF_INITIAL_SAMPLE_PROCUREMENT_OTHER: chr "" "" "" "" ...
## $ MICROMET_DETECTION_BY_IHC : chr "NO" "" "NO" "YES" ...
## $ MUTATION_COUNT : chr "24" "" "35" "18" ...
## $ NEW_TUMOR_EVENT_AFTER_INITIAL_TREATMENT : chr "" "NO" "" "" ...
## $ NTE_CENT_17_HER2_RATIO : chr "" "" "" "" ...
## $ NTE_ER_IHC_INTENSITY_SCORE : chr "" "" "" "" ...
## $ NTE_ER_STATUS : chr "" "" "" "" ...
## $ NTE_ER_STATUS_IHC_POSITIVE : chr "" "" "" "" ...
## $ NTE_HER2_FISH_STATUS : chr "" "" "" "" ...
## $ NTE_HER2_POSITIVITY_IHC_SCORE : chr "" "" "" "" ...
## $ NTE_HER2_STATUS : chr "" "" "" "" ...
## $ NTE_HER2_STATUS_IHC_POSITIVE : chr "" "" "" "" ...
## $ NTE_PR_IHC_INTENSITY_SCORE : chr "" "" "" "" ...
## $ NTE_PR_STATUS_BY_IHC : chr "" "" "" "" ...
## $ NTE_PR_STATUS_IHC_POSITIVE : chr "" "" "" "" ...
## $ OCT_EMBEDDED : chr "FALSE" "TRUE" "FALSE" "TRUE" ...
## $ ONCOTREE_CODE : chr "IMMC" "ILC" "IDC" "IDC" ...
## $ OS_MONTHS : chr "10.28" "26.02" "5.65" "40.47" ...
## $ OS_STATUS : chr "LIVING" "LIVING" "DECEASED" "LIVING" ...
## $ OTHER_PATIENT_ID : chr "DA70CF7E-0E61-4C72-B4C5-C408569D11B8" "E9002F4E-BA0B-4738-9368-DD9614A5DA8A" "d0b78f3f-a198-437a-ab8c-204345d3b75d" "161917b8-88ad-407b-ade6-d6b98478b359" ...
## $ OTHER_SAMPLE_ID : chr "025EAC6D-539E-4E43-A6F7-4A2D2DEF571D" "34E08E44-8C58-48BD-B682-1EDE7BAEAD3C" "f5675b71-ed1b-4fdb-be07-d627dc98ed88" "9ff0477c-10cd-47d6-8d0d-4e763ea185a8" ...
## $ PATHOLOGY_REPORT_FILE_NAME : chr "TCGA-A7-A3J0.A5AED5F0-A144-4615-AF78-539B8523D15F.pdf" "TCGA-OL-A66N.CDB6A529-DE8B-48F6-A569-D07251D139C3.pdf" "TCGA-AQ-A0Y5.5BC01981-366E-44EC-A986-087DB7904E5C.pdf" "TCGA-E9-A22H.672E4D81-DA47-495C-A58B-BD2B4FEA2380.pdf" ...
## $ PATHOLOGY_REPORT_UUID : chr "A5AED5F0-A144-4615-AF78-539B8523D15F" "CDB6A529-DE8B-48F6-A569-D07251D139C3" "5BC01981-366E-44EC-A986-087DB7904E5C" "672E4D81-DA47-495C-A58B-BD2B4FEA2380" ...
## $ PATH_MARGIN : chr "Negative" "Negative" "Negative" "Negative" ...
## $ PHARMACEUTICAL_TX_ADJUVANT : chr "" "NO" "" "" ...
## $ PRIMARY_SITE : chr "Right" "Left" "Right" "Left Upper Outer Quadrant" ...
## $ PROJECT_CODE : chr "" "" "" "" ...
## $ PROSPECTIVE_COLLECTION : chr "YES" "NO" "YES" "YES" ...
## $ PR_POSITIVITY_DEFINE_METHOD : chr "" "" "" "" ...
## $ PR_POSITIVITY_IHC_INTENSITY_SCORE : chr "" "" "" "" ...
## $ PR_POSITIVITY_SCALE_OTHER : chr "" "" "" "" ...
## $ PR_POSITIVITY_SCALE_USED : chr "" "" "" "" ...
## $ PR_STATUS_BY_IHC : chr "Positive" "Negative" "Positive" "Positive" ...
## $ PR_STATUS_IHC_PERCENT_POSITIVE : chr "90-99%" "" "40-49%" "10-19%" ...
## $ RACE : chr "WHITE" "WHITE" "WHITE" "WHITE" ...
## $ RADIATION_TREATMENT_ADJUVANT : chr "" "YES" "" "" ...
## $ RETROSPECTIVE_COLLECTION : chr "NO" "YES" "NO" "NO" ...
## $ SAMPLE_TYPE : chr "Primary" "Primary" "Primary" "Primary" ...
## $ SAMPLE_TYPE_ID : chr "1" "1" "1" "1" ...
## [list output truncated]Después de una vista preliminar de todos los datos que contiene nos quedamos con los que más nos interesan (explicados con detalle más adelante) y los guardamos en un dataframe.
str(df.tcga$AGE) #age## chr [1:1105] "62" "59" "70" "42" "69" "61" "50" "80" "59" "65" "79" ...str(df.tcga$MENOPAUSE_STATUS) #menostat## chr [1:1105] "Post (prior bilateral ovariectomy OR >12 mo since LMP with no prior hysterectomy)" ...str(df.tcga$AJCC_PATHOLOGIC_TUMOR_STAGE) #tgrade## chr [1:1105] "Stage IIA" "Stage IIIC" "Stage IIIA" "Stage IIB" ...str(df.tcga$LYMPH_NODE_EXAMINED_COUNT) #pnodes## chr [1:1105] "2" "15" "17" "8" "4" "4" "30" "10" "" "4" "2" "10" "12" ...str(df.tcga$ER_STATUS_BY_IHC) #strec## chr [1:1105] "Positive" "Positive" "Positive" "Positive" "Positive" ...str(df.tcga$PR_STATUS_BY_IHC) #progrec## chr [1:1105] "Positive" "Negative" "Positive" "Positive" "Positive" ...str(df.tcga$IHC_HER2) #her2## chr [1:1105] "Negative" "" "Positive" "Positive" "Equivocal" ...str(df.tcga$OS_MONTHS) #time (months)## chr [1:1105] "10.28" "26.02" "5.65" "40.47" "24.47" "97.77" "85.58" ...str(df.tcga$OS_STATUS) #cens## chr [1:1105] "LIVING" "LIVING" "DECEASED" "LIVING" "LIVING" "LIVING" ...#Remove inteterminate and empty values
#remove os_months=0
df.tcga<-df.tcga[df.tcga$OS_MONTHS>0,]
#remove menostat with values Indeterminate or empty
df.tcga<-df.tcga[as.factor(df.tcga$MENOPAUSE_STATUS)!="Indeterminate (neither Pre or Postmenopausal)",]
df.tcga<-df.tcga[as.factor(df.tcga$MENOPAUSE_STATUS)!="",]
#remove tgrade empty
df.tcga<-df.tcga[as.factor(df.tcga$AJCC_PATHOLOGIC_TUMOR_STAGE)!="",]
#remove HER, PR, ER empty or indeterminate
df.tcga<-df.tcga[as.factor(df.tcga$ER_STATUS_BY_IHC)!="",]
df.tcga<-df.tcga[as.factor(df.tcga$PR_STATUS_BY_IHC)!="",]
df.tcga<-df.tcga[as.factor(df.tcga$PR_STATUS_BY_IHC)!="Indeterminate",]
df.tcga<-df.tcga[as.factor(df.tcga$IHC_HER2)!="",]
df.tcga<-df.tcga[as.factor(df.tcga$IHC_HER2)!="Indeterminate",]
#Create dataframe
df.tcga$OS_MONTHS<-as.numeric(df.tcga$OS_MONTHS)
df.final<-data.frame(df.tcga$AGE,df.tcga$MENOPAUSE_STATUS,df.tcga$AJCC_PATHOLOGIC_TUMOR_STAGE,df.tcga$LYMPH_NODES_EXAMINED_HE_COUNT,df.tcga$ER_STATUS_BY_IHC,df.tcga$PR_STATUS_BY_IHC,df.tcga$IHC_HER2,df.tcga$OS_MONTHS,df.tcga$OS_STATUS)
#name dataframe columns
names(df.final)<-c("age", "menostat", "tgrade", "pnodes", "ER", "PR", "HER2", "time", "cens")2.2. Prepararlos para su análisis (Conversión e imputación)
Dado que muchos de los datos no están en un formato válido “”, NA,…, los imputaremos.
Para los datos numéricos, en este caso edad y número de nodos sustituiremos los valores perdidos con la media del resto.
Para los datos categóricos haremos dos secciones, aquellos que serán imputados con la moda (menostat y tgrade), y aquellos que en caso de ausencia serán interpretados como negativos (ER, HER2, PR). El valor de censura será 0 en caso de no aparecer.
impute.cols<-function(df) {
#NUMERIC->MEAN
imputar.numericos<-function(df, col) {
mean.value<-mean(as.integer(na.omit(df[,col])))
df[,col][df[,col]==""]<-floor(mean.value)
df[,col]<-as.numeric(as.character(df[,col]))
return(df[,col])
}
df$age<-imputar.numericos(df,"age")
df$pnodes<-imputar.numericos(df,"pnodes")
#CATEGORICALS
cat("LEVELS OF CATEGORICALS: \n")
#menostat
mode.menostat<-names(table(df$menostat)[which.max(table(df$menostat))])
mode.menostat
cat("Menostat levels: ", levels(df$menostat), "\n") # Peri as (Post)
df$menostat<-factor(df$menostat, labels=c("Post", "Post", "Pre"))
#tgrade
mode.tgrade<-names(table(df$tgrade)[which.max(table(df$tgrade))])
mode.tgrade
cat("Tumor grade levels: ", levels(df$tgrade), "\n") #X as mode (II)
df$tgrade<-factor(df$tgrade, labels=c("I","I", "I", "II", "II", "II", "III", "III", "III", "III", "IV", "II"))
#cens
cat("cens levels: ", levels(df.final$cens), "\n")
df$cens<-as.numeric(factor(df$cens, labels=c("1", "0"))) -1 #factor starts at 0
#ER
cat("Estrogen Receptor levels: ", levels(df$ER), "\n")
df$ER<-factor(df$ER, labels=c("Negative", "Positive"))
#HER2
cat("HER2 levels: ",levels(df$HER2), "\n") #Equivocal as Negative
df$HER2<-factor(df$HER2, labels=c("Negative", "Negative", "Positive"))
#PR
cat("Progesterone levels: ", levels(df$PR), "\n")
df$PR<-factor(df$PR, labels=c("Negative", "Positive"))
return(df)
}
df.final.imp<-impute.cols(df.final)## LEVELS OF CATEGORICALS:
## Menostat levels: Peri (6-12 months since last menstrual period) Post (prior bilateral ovariectomy OR >12 mo since LMP with no prior hysterectomy) Pre (<6 months since LMP AND no prior bilateral ovariectomy AND not on estrogen replacement)
## Tumor grade levels: Stage I Stage IA Stage IB Stage II Stage IIA Stage IIB Stage III Stage IIIA Stage IIIB Stage IIIC Stage IV Stage X
## cens levels: DECEASED LIVING
## Estrogen Receptor levels: Negative Positive
## HER2 levels: Equivocal Negative Positive
## Progesterone levels: Negative Positivedf.final.imp[1:20,]## age menostat tgrade pnodes ER PR HER2 time cens
## 1 62 Post II 0 Positive Positive Negative 10.28 1
## 2 70 Post III 5 Positive Positive Positive 5.65 0
## 3 42 Pre II 1 Positive Positive Positive 40.47 1
## 4 69 Post I 0 Positive Positive Negative 24.47 1
## 5 61 Post I 0 Positive Positive Negative 97.77 1
## 6 50 Pre III 9 Positive Positive Positive 85.58 1
## 7 80 Post II 2 Positive Positive Positive 12.19 1
## 8 59 Post III 6 Positive Positive Negative 32.98 1
## 9 65 Post II 0 Positive Positive Negative 160.78 1
## 10 49 Post II 2 Positive Negative Negative 14.78 1
## 11 62 Post II 6 Positive Negative Negative 10.58 1
## 12 64 Post II 0 Positive Positive Negative 50.82 1
## 13 64 Post II 3 Positive Positive Negative 56.90 1
## 14 58 Post II 0 Negative Positive Positive 36.17 1
## 15 85 Post II 6 Positive Positive Negative 45.04 1
## 16 79 Post II 1 Positive Positive Negative 44.78 1
## 17 54 Post III 9 Positive Positive Negative 67.81 1
## 18 36 Pre III 4 Positive Positive Positive 17.05 1
## 19 53 Post II 2 Negative Negative Negative 39.52 1
## 20 49 Pre II 0 Positive Positive Negative 13.90 1BC<-df.final.imp2.3.Datos finales
- age: Edad de la paciente en años.
#Living patients
qplot(BC$age, geom="histogram")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#Deceased patients
BC.age.nocens<-BC[BC$cens==1,]$age
qplot(BC.age.nocens, geom="histogram")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
- menostat: Estado de la menopausia, premenopausal o postmenopausal.
table(BC$menostat)##
## Post Pre
## 621 196Se encuentran en mayor proporción las pacientes Postmenopausicas.
- tgrade: Grado del tumor I, II, III, IV.
Se determina por el valor de TNM (Tumor tamaño y localización. Node ¿Se ha expandido a nodos linfáticos? ¿Donde y cuantos?. Metastasis ¿Ha habido metástasis a otras partes del cuerpo?).
Siendo:
Grado I: 20 millimeters (mm) o menor en su area más ancha.
Grado II: El tamaño del tumor es mayor a 20mm y menor a 50mm.
Grado III: Mayor a 50 mm
Grado IV: Puede ser alguno de estos casos: ha crecido en la pared del pecho(TIVa), en la piel(TIVb) o ambos(TIVc), o el cancer es inflamatorio (TIVd).
table(BC$tgrade)##
## I II III IV
## 142 484 179 12pnodes: Número de nodos linfáticos encontrados.
progrec: Receptor de progesterona (hormona esteroide C-21 involucrada en el ciclo menstrual, el embarazoy la embriogénesis) en \([fmol]\) (\(10e-15\) moles)
estrec: Receptor de estrógeno (hormonas sexuales esteroideas).
Estrógeno y Progesterona en el cancer de mama: Si las células del tumor tienen receptores de estrógeno, se dice que el cáncer es positivo en cuanto a receptores de estrógeno, sensible al estrógeno o que responde al estrógeno (ER-positivo), si no tienen, son ER-negativos y no usan el estrógeno para crecer. De forma semejante, si las células del tumor tienen receptores de progesterona, se dice que el cáncer es positivo en cuanto a receptores de progesterona (PR-positivo o PgR-positivo). Aproximadamente 80 % de los cánceres de mama tienen receptores de estrógeno (La mayoría de los cánceres de seno ER-positivos son también PR-positivos)
Si no tienen ninguno de los dos son HR negativos, receptores negativos de hormonas.
table(BC$Luminal)## < table of extent 0 >Aunque no es mucha la diferencia, podemos ver que efectivamente si hay más casos que presentan estrógeno de los que presentan progesterona.
- time: Tiempo transcurrido en estudio, en meses.
Si visualizamos el tiempo en un histograma, podremos ver claramente que no tiene una distribución (a diferencia de por ejemplo, la edad como hemos visto anteriormente), por lo que tendremos que aplicar métodos no paramétricos.
qplot(BC$time, geom="histogram")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
- cens: Si es un caso censurado(0) o no(1), es decir, si no sabemos el momento exacto de ocurrencia del evento, ya sea por abandono del estudio o porque no este no se produzca.
2.3.1 Estratificaciones de variables
2.3.1.1.Categóricas
Luminal:
Ya hemos explicado al incio de este apartado los significados de los grupos de hormonas ER, HER2 y PR.
- LuminalA: ER positivo y HER2 negativo. El grado tumoral es 1 o 2. Entre el 30-70% de los tumores de cancer de mama son de este tipo, y tienen a tener el mejor pronóstico.
-Luminal B: ER positivo. Lo tienen un 10-20% de las pacientes.
-TN (Triple negative/basal-like): ER-negative y PR-negative. Lo tienen un 15-20%, suelen ser agresivos y tener peor diagnóstico que los anteriores.
#LUMINAL
make.luminal<-function(df) {
vectorluminals <-list()
for(i in 1:nrow(df)) {
hormone_receptor_positive=(df[i,]$ER=="Positive" || df[i,]$PR=="Positive")
hormone_receptor_negative=(df[i,]$ER=="Negative" && df[i,]$PR=="Negative")
if(hormone_receptor_positive && df[i,]$HER2=="Negative") vectorluminals[i]<-"LuminalA"
else if(hormone_receptor_positive) vectorluminals[i]<-"LuminalB"
else if(df[i,]$ER=="Negative" && df[i,]$PR=="Negative" && df[i,]$HER2=="Negative") vectorluminals[i]<-"TN"
else if(hormone_receptor_negative && df[i,]$HER2=="Positive") vectorluminals[i]<-"HER2-enriched"
else if(hormone_receptor_positive && df[i,]$HER2=="Negative") vectorluminals[i]<-"Normal-line"
}
df$Luminal<-as.factor(unlist(vectorluminals))
return(df)
}
BC<-make.luminal(BC)
#NODAL AND NODE STATUS
make.nodalStatus<-function(df) {
#five groups: 0, negative; 1, one positive; 2, two to three positive; 3, four to nine positive; and 4, >9 nine positive
vectorGroups<-list()
nodeSt<-list()
for(n in 1:length(df$pnodes)) {
numnod<-df$pnodes[n]
if(numnod==0) {
vectorGroups[n]="G0"
nodeSt[n]="Node negative"
}
else if(numnod>0) {
nodeSt[n]="Node positive"
if(numnod==1) vectorGroups[n]="G1"
else if(numnod==2 || numnod==3) vectorGroups[n]="G2"
else if(numnod>=4 &&numnod<=9) vectorGroups[n]="G3"
else if(numnod>9) vectorGroups[n]="G4"
}
}
df$NodalStatus<-as.factor(unlist(vectorGroups))
df$NodeStatus<-as.factor(unlist(nodeSt))
return(df)
}
BC<-make.nodalStatus(BC)
#Otra agrupación distinta
max=max(BC$pnodes)
BC$NodalStatus2<-cut(BC$pnodes, breaks=c(0,1,3,5, max), include.lowest=T)Veamos si la estratificación de Luminal cumple las proporciones esperadas:
prop.table(table(BC$Luminal))##
## HER2-enriched LuminalA LuminalB TN
## 0.04161567 0.65850673 0.13463892 0.16523868Se adaptan perfectamente a los porcentajes descritos.
2.3.1.1.Numéricas
Estratificaremos la edad mediante clustering. También podriamos haber agrupado de esta forma el número de nodos, pero en su lugar nos hemos guiado por la definición medica con los grupos ya establecidos y lo hemos hecho de la misma manera.
Primero veremos como de bueno es el resultado para cada número de centroides, y procederemos a hacer los grupos de edad en una nueva variable.
library(factoextra)## Welcome! Related Books: `Practical Guide To Cluster Analysis in R` at https://goo.gl/13EFCZ#Data
data<-BC
data$d1<-rep(1,length(BC$age))
#Convert factor variables to numeric 'dummy' variables
data.num=model.matrix(d1+age~age, data)
#número óptimo de clusters con su puntuación silohuette
fviz_nbclust(data.frame(data$age,data$d1), kmeans, method = "silhouette") +geom_vline(xintercept = 2, linetype = 2)
Vemos que el número óptimo de clusters para nuestra variable edad según el silhouette index es 2, 6 o 10, con un valor medio de 0.6. Aún así, haremos 3 clusters ya que queremos grupos más diferenciados y también dan un resultado de clusterización bueno.
#Assign clusters
data$cl= kmeans(data.num,3)$cluster
#df<-na.omit(data.frame(data$age,data$d1))
#fviz_cluster(kmeans(data.num,3), df, geom="point") + ggtitle("k = 3")
# Vemos como se han distribuido los datos
table(data$cl)##
## 1 2 3
## 334 154 329table(data$age[data$cl==1])##
## 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70
## 19 26 12 24 26 14 26 38 26 23 19 16 14 19 17 15table(data$age[data$cl==2])##
## 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
## 15 10 11 10 9 11 9 10 14 10 6 4 4 8 5 1 2 4 3 8table(data$age[data$cl==3])##
## 26 27 28 29 30 31 32 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
## 1 1 1 4 1 1 1 8 6 8 6 9 9 11 16 10 8 9 20 22 20 21 17 25 21
## 52 53 54
## 21 22 30BC$GrupoEdadCluster<-cut(BC$age, breaks=c(25,50,66,90))
table(BC$GrupoEdadCluster)##
## (25,50] (50,66] (66,90]
## 235 363 219#Otra agrupación distinta
BC$GrupoEdad2<-cut(BC$age, breaks=c(25,37,50,90))3. Análisis
Creamos un objeto supervivencia con Surv() a partir del tiempo hasta el evento y el indicador de censura (apareceran los tiempos y un \(+\) alado de los casos que sí estén censurados).
surv.obj<-Surv(BC$time, BC$cens)
surv.obj[1:50]## [1] 10.28 5.65+ 40.47 24.47 97.77 85.58 12.19 32.98
## [9] 160.78 14.78 10.58 50.82 56.90 36.17 45.04 44.78
## [17] 67.81 17.05 39.52 13.90 31.34 54.17+ 0.33 14.72
## [25] 98.26 83.80+ 115.60 22.31 17.51 41.59 21.29 9.99
## [33] 15.44 36.83 14.72 15.47 91.92+ 16.03 30.98+ 9.00
## [41] 15.14 98.19 10.78 12.48 38.93 25.20 23.88 18.86
## [49] 44.38 82.563.1. Análisis estadístico
3.1.1. Survfit
Con survfit calculamos el estimador que se le indique, como ya hemos visto que no es una distribución paramétrica en la variable tiempo, usaremos estimadores no paramétricos.
En survfit se indica objeto surv ~ nombre de las covariables o 1, en cuyo caso lo hará para la población global, usando todas las pacientes.
Kaplain-Meier, el que usaremos principalmente es un estimador no paramétrco que calcula la supervivencia (de los pacientes cuyo evento se ha observado), en proporciones exactas.
bc.fit.km<-survfit(surv.obj~1, data=BC, type="kaplan-meier")
bc.fit.km## Call: survfit(formula = surv.obj ~ 1, data = BC, type = "kaplan-meier")
##
## n events median 0.95LCL 0.95UCL
## 817.0 738.0 27.8 24.5 31.0bc.fit.fh<-survfit(surv.obj~1, data=BC, type="fleming-harrington")
bc.fit.fh## Call: survfit(formula = surv.obj ~ 1, data = BC, type = "fleming-harrington")
##
## n events median 0.95LCL 0.95UCL
## 817.0 738.0 27.8 24.5 31.3El el objeto de tipo survfit observamos el número de pacientes, la media de supervicia y el intervalo de confianza \(95%\).
Con summary() podemos ver los datos que contiene, lo haremos para algunos tiempos en concreto.
summary(bc.fit.km, time=c(50,200))## Call: survfit(formula = surv.obj ~ 1, data = BC, type = "kaplan-meier")
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 50 222 543 0.30567 0.01665 0.274713 0.3401
## 200 1 195 0.00238 0.00235 0.000344 0.0165time el tiempo de la observación, n.risk el número de sujetos en riesgo en ese tiempo, n.evento el número de sujetos que presentaron el evento, survival la estimación de la función de supervivencia (que es menor conforme aumenta el tiempo por lo que explicaremos más adelante). Y dado que estamos estimando la población global solo con una muestra, esto conlleva un arror asociado std.err, con sus correspondientes intervalos de cofidencia CI 95%.
También tenemos la opción de tenerlo como un dataframe usando la función fortity()
bc.fit.km.df<-fortify(bc.fit.km) Aprovechamos el dataframe para observar la clara diferencia de los primeros datos con los últimos (ordenados por tiempo ascendente). Como veremos y explicaremos más adelante en las gráficas, al principio la supervivencia es prácticamente 1, así como los intervalos, y un error muy bajo, valores que cambian inversamente llegando al final de los tiempos de estudio.
bc.fit.km.df[1:10,]## time n.risk n.event n.censor surv std.err upper lower
## 1 0.03 817 1 0 0.9987760 0.001224740 1.0000000 0.9963814
## 2 0.16 816 2 0 0.9963280 0.002123916 1.0000000 0.9921891
## 3 0.23 814 1 0 0.9951040 0.002453995 0.9999018 0.9903293
## 4 0.30 813 1 0 0.9938800 0.002745339 0.9992423 0.9885466
## 5 0.33 812 15 0 0.9755202 0.005542106 0.9861744 0.9649811
## 6 0.53 797 1 0 0.9742962 0.005682535 0.9852081 0.9635051
## 7 0.62 796 1 0 0.9730722 0.005819916 0.9842354 0.9620356
## 8 0.79 795 1 0 0.9718482 0.005954464 0.9832567 0.9605722
## 9 0.99 794 4 0 0.9669523 0.006467811 0.9792880 0.9547719
## 10 1.02 790 4 0 0.9620563 0.006947980 0.9752470 0.9490440bc.fit.km.df[(nrow(bc.fit.km.df)-10):nrow(bc.fit.km.df),]## time n.risk n.event n.censor surv std.err upper
## 630 125.07 11 1 0 0.023817389 0.2705758 0.04047703
## 631 129.99 10 1 0 0.021435650 0.2903831 0.03787139
## 632 131.57 9 1 0 0.019053911 0.3133867 0.03521595
## 633 132.95 8 1 0 0.016672172 0.3406881 0.03250771
## 634 133.11 7 1 0 0.014290434 0.3740026 0.02974383
## 635 134.03 6 1 0 0.011908695 0.4161866 0.02692296
## 636 134.30 5 1 0 0.009526956 0.4724524 0.02404952
## 637 136.63 4 1 0 0.007145217 0.5536647 0.02114928
## 638 140.44 3 1 0 0.004763478 0.6879035 0.01834293
## 639 160.78 2 1 0 0.002381739 0.9865147 0.01646709
## 640 216.59 1 0 1 0.002381739 0.9865147 0.01646709
## lower
## 630 0.0140145656
## 631 0.0121328297
## 632 0.0103092926
## 633 0.0085506267
## 634 0.0068658427
## 635 0.0052675119
## 636 0.0037740005
## 637 0.0024139884
## 638 0.0012370282
## 639 0.0003444858
## 640 0.0003444858También podemos ver como sería en el caso de estratificarse con alguna variable, menostat por ejemplo, la que mayor diferencia muestra a priori.
bc.fit.km.menostat<-survfit(surv.obj~menostat, data=BC, type="kaplan-meier")
head(fortify(bc.fit.km.menostat, strata=menostat))## time n.risk n.event n.censor surv std.err upper lower
## 1 0.03 621 1 0 0.9983897 0.001611604 1.0000000 0.9952411
## 2 0.16 620 2 0 0.9951691 0.002795893 1.0000000 0.9897306
## 3 0.30 618 1 0 0.9935588 0.003231035 0.9998707 0.9872867
## 4 0.33 617 13 0 0.9726248 0.006732249 0.9855436 0.9598753
## 5 0.53 604 1 0 0.9710145 0.006933169 0.9842994 0.9579089
## 6 0.62 603 1 0 0.9694042 0.007129068 0.9830445 0.9559532
## strata
## 1 Post
## 2 Post
## 3 Post
## 4 Post
## 5 Post
## 6 Postggsurvplot(bc.fit.km, xlab = "Tiempo", censor = T, conf.int=T, ylab = "Survival Probability", title = "Survival probability", ggtheme = theme_bw())
ggsurvplot(bc.fit.km, xlab = "Tiempo", censor = T, conf.int=T, ylab = "Survival Probability", title = "Survival probability", ggtheme = theme_bw(), surv.median.line = "h", risk.table = T, cumevents=T)
Junto con la curva, también tenemos la opción de representar las tablas de riesgo y eventos de esta misma.
ggsurvplot(bc.fit.km.menostat, xlab = "Tiempo", censor = T, conf.int=T, ylab = "Survival Probability", title = "Survival probability", pval=TRUE, ggtheme = theme_bw())
3.1.2.Desviación estandar e intervalos de confianza
Se definen directamente en la función survfit().
Se puede definir el error: tipo de estimación para las desviaciones, conf.type para definit el tipo de transformación para calcular los intervalos de confianza y conf.int, el nivel de confianza (95% por defecto) (Este nivel quiere decir que ese tanto porciento de las veces la media se encontrará dentro del intervalo).
error: greenwood(defecto), tsiatisconf.type: plain, log(defecto), log-log
surv.tsiatis.plain<-survfit(surv.obj ~ 1, data=BC, error="tsiatis", conf.type = "plain", conf.int = 0.99)
head(surv.tsiatis.plain)## Call: survfit(formula = surv.obj ~ 1, data = BC, error = "tsiatis",
## conf.type = "plain", conf.int = 0.99)
##
## n events median 0.99LCL 0.99UCL
## 817.0 738.0 27.8 23.7 32.3surv.tsiatis.loglog<-survfit(surv.obj ~ 1, data=BC, error="tsiatis", conf.type = "log-log", conf.int = 0.90)
head(surv.tsiatis.loglog)## Call: survfit(formula = surv.obj ~ 1, data = BC, error = "tsiatis",
## conf.type = "log-log", conf.int = 0.9)
##
## n events median 0.9LCL 0.9UCL
## 817.0 738.0 27.8 24.7 30.63.1.3.Media, mediana, percentiles
media = \(\int^{\tau}_{0}S(t)dt\) siendo S(t) la función de supervivencia estimada hasta valor \(\tau\).
print(bc.fit.km, print.rmean=T)## Call: survfit(formula = surv.obj ~ 1, data = BC, type = "kaplan-meier")
##
## n events *rmean *se(rmean) median 0.95LCL
## 817.00 738.00 40.62 1.27 27.83 24.47
## 0.95UCL
## 30.98
## * restricted mean with upper limit = 217probs<-c(0.05, 0.5, 0.95)
quantile(bc.fit.km, probs)## $quantile
## 5 50 95
## 2.60 27.83 108.28
##
## $lower
## 5 50 95
## 1.02 24.47 105.98
##
## $upper
## 5 50 95
## 5.85 30.98 118.50quantile(bc.fit.km.menostat, probs)## $quantile
## 5 50 95
## menostat=Post 1.77 24.11 107.85
## menostat=Pre 7.29 40.37 115.18
##
## $lower
## 5 50 95
## menostat=Post 0.99 22.08 101.54
## menostat=Pre 4.40 32.75 107.62
##
## $upper
## 5 50 95
## menostat=Post 5.32 27.56 118.36
## menostat=Pre 11.33 50.33 NA3.1.4. Estimación de la función de riesgo acumulado
Nos permitirá ver como aumenta con el tiempo el riesgo de que ocurra el evento. Como es de esperar, veremos claramente como es inversa a la curva de supervivencia, ya que cuando el riesgo sube, esta baja.
Tenemos varias formas de estimarla, entre ellas exploraremos máxima verosimilitud y Nelson-Aalen
El riesgo acumulado vendrá dado por CumHaz= cumsum(n.event/n.risk), usando la función mutate (perteneciente a la libreria dplyr).Para visualizarlo simplemente tendremos que indicar fun="cumhaz".
Podemos visualizarlo para los diferentes niveles de un evento en concreto usando group_by.
3.1.4.1. Máxima verosimilitud
En este caso se estima igual que anteriormente, pero se calcula la suma acumulada del número de enventos entre las personas en riesgo.
bc.fit.ra<-survfit(surv.obj~1) %>% fortify %>% mutate(CumHaz = cumsum(n.event/n.risk))
head(bc.fit.ra)## time n.risk n.event n.censor surv std.err upper lower
## 1 0.03 817 1 0 0.9987760 0.001224740 1.0000000 0.9963814
## 2 0.16 816 2 0 0.9963280 0.002123916 1.0000000 0.9921891
## 3 0.23 814 1 0 0.9951040 0.002453995 0.9999018 0.9903293
## 4 0.30 813 1 0 0.9938800 0.002745339 0.9992423 0.9885466
## 5 0.33 812 15 0 0.9755202 0.005542106 0.9861744 0.9649811
## 6 0.53 797 1 0 0.9742962 0.005682535 0.9852081 0.9635051
## CumHaz
## 1 0.001223990
## 2 0.003674971
## 3 0.004903472
## 4 0.006133484
## 5 0.024606391
## 6 0.025861096ggsurvplot(bc.fit.ra, fun="cumhaz", xlab="Tiempo(meses)", ylab="Riempo acumulado", main="Riesgo acumulado", censor=T, conf.int=T, palette="Dark2", ggtheme = theme_bw())
bc.fit.ra.menostat<-survfit(surv.obj~menostat, data=BC) %>% fortify %>% group_by(strata) %>% mutate(CumHaz = cumsum(n.event/n.risk))
head(bc.fit.ra.menostat)## # A tibble: 6 x 10
## # Groups: strata [1]
## time n.risk n.event n.censor surv std.err upper lower strata CumHaz
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <dbl>
## 1 0.03 621 1 0 0.998 0.00161 1 0.995 Post 0.00161
## 2 0.16 620 2 0 0.995 0.00280 1 0.990 Post 0.00484
## 3 0.3 618 1 0 0.994 0.00323 1.000 0.987 Post 0.00645
## 4 0.33 617 13 0 0.973 0.00673 0.986 0.960 Post 0.0275
## 5 0.53 604 1 0 0.971 0.00693 0.984 0.958 Post 0.0292
## 6 0.62 603 1 0 0.969 0.00713 0.983 0.956 Post 0.0308ggsurvplot(survfit(surv.obj~menostat, data=BC), fun="cumhaz", xlab="Tiempo(meses)", ylab="Riesgo acumulado", main="Riesgo acumulado", censor=T, conf.int=T, legend.title="Perfil menostat", palette="Dark2", ggtheme = theme_bw())
También podemos ver los valores de riesgo en una tabla:
ggsurvplot(bc.fit.km.menostat, xlab = "Tiempo", censor = T, conf.int=T, ylab = "Survival Probability", title = "Survival probability", pval=TRUE, risk.table = TRUE, risk.table.y.text.col = TRUE, ggtheme = theme_bw())
3.1.4.2. Nelson-Aalaen
La información es extraída de la misma manera que la Estimación Máxima Verosimilitud, con la diferencia que se calcula la suma acumulada del número de eventos entre las personas en riesgo.
qplot(time, CumHaz, data=bc.fit.ra, geom="step", xlab="Tiempo(meses)", ylab="Riesgo acumulado", main="Riesgo acumulado")
qplot(time, CumHaz, col=strata, data=bc.fit.ra.menostat, geom="step", xlab="Tiempo(meses)", ylab="Riesgo acumulado", main="Riesgo acumulado")
3.2. Métodos no paramétricos
Todo lo visto anteriormente se realizada de forma no paramétrica, pero solo lo hemos visto para Kaplain-Meier.
Veamos si otros métodos nos dan alguna diferencia.
km<-survfit(surv.obj~1, data=BC, type="kaplan-meier")
fh<-survfit(surv.obj~1, data=BC, type="fleming-harrington")
fh2<-survfit(surv.obj~1, data=BC, type="fh2")
list.fit=list(km,fh)
ggplot() + geom_step(aes(time, surv, col = "Kaplain-meier"), data = km) + geom_step(aes(time, surv, col = "Fleming-harrington"), data = fh) + geom_step(aes(time, surv, col = "fh2"), data=fh2) + labs(x = "Tiempo (meses)", y = "Probabilidad de Supervivencia", col = "Ajustes no parametricos", title = "Comparacion de ajustes en curvas de supervivencia")
La diferencia de Kaplain-meier con otros ajustes es inapreciable a simple vista, por lo que seguiremos usando este primero.
3.3. Métodos paramétricos
Estas técnicas a diferencia de las no paramétricas, tienen que ajustarse a una distribución de datos, y cumplir condiciones de validez.
Para estimar las curvas de supervivencia de forma no paramétrica, usaremos la función flexsurvreg(), donde los parámetros se estiman por máxima verosimilitud (visto en 3.1.4.1).
Tenemos disponibles multitud de distribuciones por las que ajustar (indicadas en el parámetro dist), algunas de ellas son:
| Nombre | dist | Numero de parámetros |
|---|---|---|
| Gamma generalizada (stable) | “gengamma” | 3 |
| Gamma generalizada (original) | “gengamma.orig” | 3 |
| F generalizada (stable) | “genf” | 4 |
| F generalizada (original) | “genf.orig” | 4 |
| Weibull | “weibull” | 2 |
| Gamma | “gamma” | 2 |
| Exponencial | “exp” | 1 |
| Log-logistica | “llogis” | 2 |
| Log-normal | “lnorm” | 2 |
| Gompertz | “gompertz” | 2 |
Los parámetros definidos como positivos se estiman con la logarítmica y los intervalos de confianza se calculan a partir de la matriz Hessiana.
Compararemos graficamente algunas de ellas aplicadas a nuestro dataset:
#Distribuciones
dist<-c("exp", "gengamma", "gamma", "weibull", "llogis", "lnorm", "gompertz", "genf")
#Lista de ajustes de las distribuciones
fit.param<-sapply(dist, function(x) flexsurvreg(surv.obj~1, dist=x), USE.NAMES=T, simplify=F)
param.dist<-names(fit.param) %>% sapply(function(x) cbind(dist = x, data.frame(summary(fit.param[[x]]))), simplify = F) %>% do.call("rbind", .)
x <- (1:length(dist)) + 1 #Colores
names(x) <- dist
#Supervivence
ggplot() + geom_line(aes(time, est, col = factor(dist)), param.dist, size=0.9, alpha=0.5) + labs(col = "Distribuciones parametricas", linetype = "Ajustes", x ="Tiempo(meses)", y = "Probabilidad de Supervivencia", title = "Comparacion de ajustes parametricos")
#Risk
ggplot() + geom_line(aes(time, -log(est), col = factor(dist)), param.dist, size=0.9, alpha=0.5) + labs(col = "Distribuciones parametricas", linetype = "Ajustes", x ="Tiempo(meses)", y = "Riesgo acumulado", title = "Comparacion de ajustes parametricos")
También podemos ver los valores de AIC, BIC y LogLik para cada distribución.
AIC (Akaike information criterion): Estimador de la calidad relativa del modelo estadístico para el dataset al que se ha aplicado. Representa la cantidad de información perdida, por tanto el modelo será mejor para un menor AIC.
BIC (Bayesian information criterion): Criterio para seleccionar alguno de los modelos, es preferible el de menos BIC.
LogLik (likelihood function)
AIC.model <- sapply(fit.param, function(x) c(AIC = AIC(x), BIC = BIC(x), LogLik = logLik(x)), simplify = T)
AIC.model## exp gengamma gamma weibull llogis lnorm
## AIC 6972.944 6957.689 6958.970 6956.217 7027.395 7155.689
## BIC 6977.649 6971.806 6968.382 6965.629 7036.806 7165.101
## LogLik -3485.472 -3475.844 -3477.485 -3476.109 -3511.697 -3575.845
## gompertz genf
## AIC 6959.958 6959.691
## BIC 6969.369 6978.513
## LogLik -3477.979 -3475.845No hay demasiada diferencia entre los modelos, pero a simple vista el de weibull parace ser el más adecuado en este caso.
3.4.Regresión de Cox
El objetivo principal de Cox (modelo semiparamétrico de regresión multivariante) es poder estimar el efecto de riesgo conjunto de distintas variables, además, funciona con variables continuas además de categóricas, permitiendonos esto analizar otras variables no vistas anteriormente como por ejemplo la edad.
Podemos ver como cada factor influye al Hazard rate o failure rate. Como cambia la supervivencia respecto a la otra variable. Es la probabilidad de que si algo sobrevive a un momento, lo hará también para el siguiente.
Función de riesgo: \(h(t)=\frac{f(t)}{R(t)}\) (riesgo de que se produzca en evento en el instante de tiempo \(t\).)
Siendo:
\(f(t)\): Función de probabilidad de densidad (probabilidad de que el evento caiga en un intervalo específico).
\(R(t)\): Función de supervivencia dado por las covariantes \(x_1, x_2,... x_p\) indicadas, siendo estos multiplicados por sus respectivos impactos \(b_1, b_2,...b_p\).
\(h0\): Baseline hazard (no covariables, all \(x_i\)=0)
Hazard ratios (HR) \(exp(b_i)\)
- HR=1: Esa variable no tiene efecto
- HR <1: Reducción de la probabilidad hazard
- HR >1: aumento
Prueba de hazards proporcionales o proporcionalidad de los riesgos:
“If an individual has a risk of death at some initial time point that is twice as high as that of another individual, then at all later times the risk of death remains twice as high.” Se realizarán con cox.zph para aceptar el modelo como válido, si esta da un \(p-value \ge 0.05\). De esta forma sabemos que se mantiene la diferencia entre las pacientes.
3.4.1. Busqueda exhaustiva
3.4.1.1. Modelos univariantes
library(FSelector)
returnSignificatives<-function(vbles) {
variables<-list()
pvalues<-list()
for (v in vbles) {
form<-as.simple.formula(v,"surv.obj")
sum<-summary(coxph(form, data=BC))$coefficients
p<-sum[,"Pr(>|z|)"]
if(p<0.05){
variables=c(variables,v)
pvalues=c(pvalues,p)
}
}
significancia<-data.frame(unlist(variables),unlist(pvalues))
names(significancia)<-c("variable", "p-value")
return(significancia)
}
returnSignificatives(c("age", "menostat","tgrade","pnodes", "Luminal", "NodalStatus", "NodeStatus"))## variable p-value
## 1 age 7.483853e-05
## 2 menostat 1.362899e-04
## 3 pnodes 2.374573e-03Parece ser que individualmente, las más significativas son edad, menopausia y número de nodos.
3.4.1.2. Modelos multivariantes
Probaremos juntas las más significativas, para ver si explican mejor la función de riesgo, pero algunas como NodeStatus o tgrade no las descartaremos aún ya que pueden dar buenos resultados.
Al añadir más variables, se tomará una de ellas como referencia.
Con la función de búsqueda exhaustiva intentaremos encontrar la mejor combinación posible de variables que expliquen la función de riesgo. Nos quedaremos en cada nivel con la mejor (la de menor \(p-value\)) e iremos añadiendo otras solo en el caso de que su incorporación mejore el resultado que ya teníamos. Nos quedamos con la evolución de los \(p-values\).
exhaustiveSearch <-function(lista.vbles) {
pvalue.acumulado<-list()
for(nivel in 1:length(lista.vbles)) { #cada nivel (#vbles=#nivel)
best.pvalue<-1
for(v in lista.vbles) { #vbles de cada nivel
ifelse(nivel==1, att<-v, att<-c(v, lista.acumulada))
ec <-as.simple.formula(att,"surv.obj")
sum<-summary(coxph(ec, data=BC))
new.p<-sum$logtest["pvalue"]
if(new.p <0.05 && new.p<best.pvalue) {
best.pvalue<-new.p
best.vble.nivel<-v
}
}
lista.vbles<-lista.vbles[lista.vbles != best.vble.nivel]
ifelse(nivel!=1, pvalue.acumulado<-c(pvalue.acumulado, best.pvalue), first.pvalue<-best.pvalue)
ifelse(nivel==1, lista.acumulada<-best.vble.nivel,lista.acumulada<-c(lista.acumulada, best.vble.nivel))
ifelse(nivel==1, p.acum.ac<-best.pvalue, p.acum.ac<-pvalue.acumulado[nivel-1])
if(length(lista.acumulada)!=nivel || p.acum.ac>best.pvalue[1] || length(lista.acumulada)==10) break
}
df<-data.frame(lista.acumulada,unlist(list(c(first.pvalue,pvalue.acumulado))))
names(df)<-c("vble acumulada", "pvalue_acumulado")
return(df)
}
es<-exhaustiveSearch(c("age", "menostat","tgrade","pnodes", "Luminal", "NodalStatus", "NodalStatus2", "NodeStatus", "GrupoEdadCluster", "GrupoEdad2"))
df.order<-es[order(es$pvalue_acumulado),]
df.order## vble acumulada pvalue_acumulado
## 3 tgrade 5.956164e-09
## 4 pnodes 9.276238e-09
## 5 GrupoEdad2 1.690446e-08
## 2 menostat 3.034877e-08
## 6 age 3.851068e-08
## 7 GrupoEdadCluster 8.173465e-08
## 8 NodeStatus 1.851588e-07
## 9 NodalStatus 3.260355e-07
## 10 Luminal 2.267670e-06
## 1 NodalStatus2 6.423682e-06Las variables más influyentes son tgrade, pnodes y GrupoEdad. Efectivamente coincide con las variables encontrados en el análisis multivariante, pero ahora el pvalue es mucho menor, pasando de ser del orden \(e-05, e-04, e-03\) a \(e-09\), lo cual nos confirma que juntas explican mejor la función de rieso.
En el primer puesto tenemos a tgrade, variable que no aparecia anteriormente. Sola no es significativa, pero sí en conjunto.
Así como para añadir una variable a analizar se usa +, para ver todas las interacciones posibles entre las variables usaremos *.
Comprobamos que sea válido antes de examinarlo.
El coxph() nos dará la siguiente información:
- coef: Coeficiente estimado de la beta (los valores beta, en función del valor explican la importancia de esta variable para la supervivencia)
- exp(coef): Exponencial elevado al coeficiente beta estimado.
- se(coef): La desviación estándar de la estimación.
- z: Valor del estadístico para prueba Wald de beta igual a cero.
- p: p-valor de la prueba Wald, beta igual a cero.
Usando summary() los resultados son una versión más detallada de la función coxph, en esta se encuentran agregada una tabla con los coeficientes estimados de los exponentes β con sus respectivos intervalos de confianza.
Además se agregan estadísticos de Concordancia y R2; así como los estadísticos de prueba, grados de libertad y p-valores de las pruebas de hipótesis de Razón de Verosimilitud, Prueba de Wald y Prueba de Puntajes (score test).
res.cox <- coxph(surv.obj~tgrade+pnodes+GrupoEdadCluster, data=BC)
cox.zph(res.cox, transform="km", global=TRUE) #ACEPTADO## rho chisq p
## tgradeII -0.10345 7.8452 0.0051
## tgradeIII -0.06237 2.9899 0.0838
## tgradeIV -0.05267 2.0735 0.1499
## pnodes 0.00662 0.0288 0.8653
## GrupoEdadCluster(50,66] -0.03672 0.9831 0.3214
## GrupoEdadCluster(66,90] -0.07404 4.0399 0.0444
## GLOBAL NA 12.1104 0.0596summ<-summary(res.cox)
summ## Call:
## coxph(formula = surv.obj ~ tgrade + pnodes + GrupoEdadCluster,
## data = BC)
##
## n= 817, number of events= 738
##
## coef exp(coef) se(coef) z Pr(>|z|)
## tgradeII 0.09981 1.10496 0.09961 1.002 0.316336
## tgradeIII 0.05582 1.05741 0.14269 0.391 0.695658
## tgradeIV -0.98886 0.37200 0.51252 -1.929 0.053678 .
## pnodes 0.02798 1.02837 0.01089 2.570 0.010182 *
## GrupoEdadCluster(50,66] 0.28066 1.32401 0.08752 3.207 0.001342 **
## GrupoEdadCluster(66,90] 0.37929 1.46125 0.10241 3.703 0.000213 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## tgradeII 1.105 0.9050 0.9090 1.343
## tgradeIII 1.057 0.9457 0.7994 1.399
## tgradeIV 0.372 2.6882 0.1362 1.016
## pnodes 1.028 0.9724 1.0067 1.051
## GrupoEdadCluster(50,66] 1.324 0.7553 1.1153 1.572
## GrupoEdadCluster(66,90] 1.461 0.6843 1.1955 1.786
##
## Concordance= 0.578 (se = 0.012 )
## Rsquare= 0.037 (max possible= 1 )
## Likelihood ratio test= 31.06 on 6 df, p=2e-05
## Wald test = 30.73 on 6 df, p=3e-05
## Score (logrank) test = 31.13 on 6 df, p=2e-05res.cox1 <- coxph(surv.obj~tgrade+NodeStatus+NodalStatus, data=BC)
cox.zph(res.cox1, transform="km", global=TRUE) #DESCARTADO## rho chisq p
## tgradeII -0.091220 6.01e+00 0.0142
## tgradeIII -0.025336 5.07e-01 0.4767
## tgradeIV -0.039975 1.22e+00 0.2702
## NodeStatusNode positive 0.020980 3.45e-01 0.5571
## NodalStatusG1 -0.023330 4.32e-01 0.5112
## NodalStatusG2 -0.000196 2.98e-05 0.9956
## NodalStatusG3 -0.084237 5.23e+00 0.0222
## NodalStatusG4 NA NaN NaN
## GLOBAL NA 2.35e+01 0.0028res.cox2 <- coxph(surv.obj~menostat+tgrade+NodeStatus+Luminal, data=BC)
cox.zph(res.cox2, transform="km", global=TRUE) #DESCARTADO## rho chisq p
## menostatPre 0.094968 6.58e+00 0.0103
## tgradeII -0.080805 4.77e+00 0.0289
## tgradeIII -0.043925 1.43e+00 0.2314
## tgradeIV -0.049536 1.84e+00 0.1755
## NodeStatusNode positive -0.045546 1.55e+00 0.2130
## LuminalLuminalA 0.010433 8.15e-02 0.7753
## LuminalLuminalB 0.026295 5.13e-01 0.4737
## LuminalTN -0.000543 2.19e-04 0.9882
## GLOBAL NA 1.70e+01 0.0300res.cox3 <- coxph(surv.obj~menostat+tgrade+NodeStatus+Luminal+GrupoEdadCluster, data=BC)
cox.zph(res.cox3, transform="km", global=TRUE) #ACEPTADO## rho chisq p
## menostatPre 0.07173 3.4189 0.0645
## tgradeII -0.08422 5.2189 0.0223
## tgradeIII -0.04750 1.6938 0.1931
## tgradeIV -0.05252 2.0714 0.1501
## NodeStatusNode positive -0.04050 1.2398 0.2655
## LuminalLuminalA 0.00593 0.0267 0.8701
## LuminalLuminalB 0.02757 0.5663 0.4517
## LuminalTN -0.00431 0.0139 0.9062
## GrupoEdadCluster(50,66] 0.01945 0.2555 0.6132
## GrupoEdadCluster(66,90] -0.01017 0.0709 0.7900
## GLOBAL NA 17.9113 0.0565summ3<-summary(res.cox3)
summ3## Call:
## coxph(formula = surv.obj ~ menostat + tgrade + NodeStatus + Luminal +
## GrupoEdadCluster, data = BC)
##
## n= 817, number of events= 738
##
## coef exp(coef) se(coef) z Pr(>|z|)
## menostatPre -0.213294 0.807918 0.122083 -1.747 0.0806 .
## tgradeII 0.059617 1.061430 0.106317 0.561 0.5750
## tgradeIII 0.137181 1.147036 0.141452 0.970 0.3321
## tgradeIV -0.959077 0.383247 0.513807 -1.867 0.0620 .
## NodeStatusNode positive 0.125685 1.133925 0.089555 1.403 0.1605
## LuminalLuminalA -0.149212 0.861387 0.193617 -0.771 0.4409
## LuminalLuminalB -0.004047 0.995961 0.212315 -0.019 0.9848
## LuminalTN -0.120004 0.886917 0.209885 -0.572 0.5675
## GrupoEdadCluster(50,66] 0.171218 1.186750 0.118177 1.449 0.1474
## GrupoEdadCluster(66,90] 0.246775 1.279892 0.134044 1.841 0.0656 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## menostatPre 0.8079 1.2377 0.6360 1.026
## tgradeII 1.0614 0.9421 0.8618 1.307
## tgradeIII 1.1470 0.8718 0.8693 1.513
## tgradeIV 0.3832 2.6093 0.1400 1.049
## NodeStatusNode positive 1.1339 0.8819 0.9514 1.351
## LuminalLuminalA 0.8614 1.1609 0.5894 1.259
## LuminalLuminalB 0.9960 1.0041 0.6569 1.510
## LuminalTN 0.8869 1.1275 0.5878 1.338
## GrupoEdadCluster(50,66] 1.1867 0.8426 0.9414 1.496
## GrupoEdadCluster(66,90] 1.2799 0.7813 0.9842 1.664
##
## Concordance= 0.577 (se = 0.012 )
## Rsquare= 0.039 (max possible= 1 )
## Likelihood ratio test= 32.64 on 10 df, p=3e-04
## Wald test = 30.93 on 10 df, p=6e-04
## Score (logrank) test = 31.32 on 10 df, p=5e-04#cox1
summ$logtest[3]## pvalue
## 2.468529e-05#cox3
summ3$logtest[3]## pvalue
## 0.0003134935res.cox4 <- coxph(surv.obj~tgrade+NodeStatus+Luminal+GrupoEdadCluster, data=BC)
cox.zph(res.cox4, transform="km", global=TRUE) #ACEPTADO## rho chisq p
## tgradeII -0.08405 5.1935 0.0227
## tgradeIII -0.04689 1.6479 0.1992
## tgradeIV -0.05256 2.0742 0.1498
## NodeStatusNode positive -0.03988 1.2030 0.2727
## LuminalLuminalA 0.00814 0.0505 0.8223
## LuminalLuminalB 0.02946 0.6455 0.4217
## LuminalTN -0.00422 0.0133 0.9082
## GrupoEdadCluster(50,66] -0.03870 1.1224 0.2894
## GrupoEdadCluster(66,90] -0.07259 3.8113 0.0509
## GLOBAL NA 15.0513 0.0895summ4<-summary(res.cox4)
summ4$logtest[3]## pvalue
## 0.0005183621| Variables | p-value propircionalidad | p-value modelo |
|---|---|---|
| tgrade+pnodes+GrupoEdadCluster | 0.0596 | 2.468529e-05 |
| menostat+tgrade+NodeStatus+Luminal+GrupoEdadCluster | 0.0565 | 0.0003134935 |
| tgrade+NodeStatus+Luminal+GrupoEdadCluster | 0.0895 | 0.0005183621 |
Algunos daban un buen p-value como modelo pero hemos tenido que descartarlo ya que no cumplia la proporcionalidad. De entre los restantes probados, el segundo es el más fiable en cuanto a proporcionalidad, pero con poca diferencia respecto al primero, el cual marca una notable diferencia en cuanto al p-value del modelo de cox.
El resultado de proporcionalidad mejora mucho al quitar menostat, pero no el pvalue del modelo.
Probaremos alguno más siguiendo el patrón con mejores resultados.
res.cox5 <- coxph(surv.obj~menostat+tgrade+NodeStatus+NodalStatus+Luminal+GrupoEdadCluster, data=BC)
cox.zph(res.cox5, transform="km", global=TRUE) #RECHAZADO## rho chisq p
## menostatPre 0.06805 2.97541 0.08454
## tgradeII -0.09666 6.85358 0.00885
## tgradeIII -0.03388 0.93281 0.33413
## tgradeIV -0.04972 1.89922 0.16817
## NodeStatusNode positive 0.02557 0.52991 0.46664
## NodalStatusG1 -0.03067 0.76689 0.38118
## NodalStatusG2 -0.00489 0.01898 0.89041
## NodalStatusG3 -0.08658 5.60686 0.01789
## NodalStatusG4 NA NaN NaN
## LuminalLuminalA 0.00139 0.00147 0.96940
## LuminalLuminalB 0.03365 0.84761 0.35723
## LuminalTN -0.00643 0.03108 0.86007
## GrupoEdadCluster(50,66] 0.01962 0.25098 0.61639
## GrupoEdadCluster(66,90] -0.01056 0.07693 0.78150
## GLOBAL NA 31.41162 0.00485Al añador el grupo nodal se anula el criterio de proporcionalidad.
res.cox6 <- coxph(surv.obj~menostat+NodeStatus+Luminal+GrupoEdadCluster, data=BC)
cox.zph(res.cox6, transform="km", global=TRUE) #ACEPTADO## rho chisq p
## menostatPre 0.07539 3.78033 0.0519
## NodeStatusNode positive -0.07990 4.74086 0.0295
## LuminalLuminalA 0.01216 0.11216 0.7377
## LuminalLuminalB 0.02473 0.45645 0.4993
## LuminalTN -0.00315 0.00743 0.9313
## GrupoEdadCluster(50,66] 0.01726 0.20065 0.6542
## GrupoEdadCluster(66,90] -0.00680 0.03171 0.8587
## GLOBAL NA 12.09198 0.0976summ6<-summary(res.cox6)Al quitar el grado mejora mucho la proporcionalidad.
survfit(surv.obj ~ tgrade, BC, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia por grado", conf.int = T, legend.title = "Grado", pval=TRUE, ggtheme = theme_bw(), risk.table=T, risk.table.col="strata", risk.table.height = 0.32)
Si visualizamos su curva de supervivencia podemos ver como tiene un pvalue poco fiable.El nodo IV tiene la menor supervivencia como podemos ver muy claramente, pero al haber tan pocos pacientes en ese grupo sus intervalos de confianza son muy amplios para poder alcanzar el 95% que le hemos marcado por defecto, solapando completamente al resto de niveles.
res.cox7 <- coxph(surv.obj~menostat+NodeStatus+NodalStatus+Luminal+GrupoEdadCluster, data=BC)
cox.zph(res.cox7, transform="km", global=TRUE) #ACEPTADO## Warning in cor(xx, r2): the standard deviation is zero## rho chisq p
## menostatPre 0.07219 3.4030 0.0651
## NodeStatusNode positive 0.02766 0.5723 0.4493
## NodalStatusG1 -0.04701 1.6636 0.1971
## NodalStatusG2 -0.01631 0.1994 0.6552
## NodalStatusG3 -0.09283 6.4465 0.0111
## NodalStatusG4 NA NaN NaN
## LuminalLuminalA 0.00649 0.0321 0.8578
## LuminalLuminalB 0.02755 0.5680 0.4511
## LuminalTN -0.00547 0.0225 0.8809
## GrupoEdadCluster(50,66] 0.01659 0.1826 0.6691
## GrupoEdadCluster(66,90] -0.00624 0.0268 0.8699
## GLOBAL NA 23.9073 0.0131summ7<-summary(res.cox7)
res.cox8 <- coxph(surv.obj~NodeStatus+NodalStatus+Luminal, data=BC)
summ8<-summary(res.cox8)Habiendo quitado el grado tumoral que causaba problemas, se puede añadir NodalStatus y la propocionalidad mejora bastante.
summ4$logtest[3]## pvalue
## 0.0005183621summ6$logtest[3]## pvalue
## 0.0006044408summ7$logtest[3]## pvalue
## 0.0002492881summ8$logtest[3]## pvalue
## 0.0430766A diferencia de los demás que no marcan gran diferencia, quitarle la información sobre el grado tumoral el resultado mejora enormemente, además, la edad y estado menopausico mejoran mucho el pvalue, por lo que parecen ser importantes en el estudio. Al añadirle otros empeora un poco el resultado por lo que en principio menostat+NodeStatus+NodalStatus+Luminal+GrupoEdadCluster serán neustras covariantes principales.
Exploramos varios valores y finalmente podemos afinar los levels que mejor explican la función de riesgo, dando el menor p-value encontrado ahasta el momento:
res.cox8<-coxph(surv.obj~+I(menostat=="Post")+ I(NodeStatus=="NodePositive") + I(NodalStatus=="G2")+I(GrupoEdadCluster=="(50,66]")+I(Luminal=="TN"), data=BC)
summ8<-summary(res.cox8)
summ8$logtest[3]## pvalue
## 0.003594444summ7## Call:
## coxph(formula = surv.obj ~ menostat + NodeStatus + NodalStatus +
## Luminal + GrupoEdadCluster, data = BC)
##
## n= 817, number of events= 738
##
## coef exp(coef) se(coef) z Pr(>|z|)
## menostatPre -0.22285 0.80023 0.12101 -1.842 0.0655 .
## NodeStatusNode positive 0.31993 1.37703 0.18799 1.702 0.0888 .
## NodalStatusG1 -0.31326 0.73106 0.20566 -1.523 0.1277
## NodalStatusG2 -0.26269 0.76898 0.20506 -1.281 0.2002
## NodalStatusG3 -0.01144 0.98863 0.19368 -0.059 0.9529
## NodalStatusG4 NA NA 0.00000 NA NA
## LuminalLuminalA -0.13895 0.87027 0.19366 -0.717 0.4731
## LuminalLuminalB -0.04109 0.95974 0.21239 -0.193 0.8466
## LuminalTN -0.10574 0.89966 0.21011 -0.503 0.6148
## GrupoEdadCluster(50,66] 0.13973 1.14996 0.11722 1.192 0.2332
## GrupoEdadCluster(66,90] 0.19479 1.21505 0.13409 1.453 0.1463
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## menostatPre 0.8002 1.2496 0.6313 1.014
## NodeStatusNode positive 1.3770 0.7262 0.9526 1.990
## NodalStatusG1 0.7311 1.3679 0.4885 1.094
## NodalStatusG2 0.7690 1.3004 0.5145 1.149
## NodalStatusG3 0.9886 1.0115 0.6763 1.445
## NodalStatusG4 NA NA NA NA
## LuminalLuminalA 0.8703 1.1491 0.5954 1.272
## LuminalLuminalB 0.9597 1.0419 0.6329 1.455
## LuminalTN 0.8997 1.1115 0.5960 1.358
## GrupoEdadCluster(50,66] 1.1500 0.8696 0.9139 1.447
## GrupoEdadCluster(66,90] 1.2151 0.8230 0.9342 1.580
##
## Concordance= 0.584 (se = 0.012 )
## Rsquare= 0.04 (max possible= 1 )
## Likelihood ratio test= 33.23 on 10 df, p=2e-04
## Wald test = 33.3 on 10 df, p=2e-04
## Score (logrank) test = 33.59 on 10 df, p=2e-04Los Luminals y Grupos nodales tienen el valor ß negativo , lo cual quiere decir que a mayor sea este, menor es el riesgo, los grupos de edad por el contrario tiene un ß positivo, a mayor número de nodos mayor riesgo, y de esta forma interpretamos el riesgo.
Los evaluaremos más adelante, pero podemos intuir que el riesgo principalmente para el estado nodal positivo.
También volvemos a ver que la (pre)menopausia influye enormemente, llegando a ser la de menor p-value.
El intervalo de confianza con mayor rango lo tenemos en el estado nodal.
Otro efecto que también hemos observado es que a mayor edad, menor es el riesgo.
3.4.2. Comparación modelos
La comparación entre los modelos se realiza con la función anova(), quw realiza la prueba de Razón de Verosimilitud entre dos modelo, le daremos como argumentos de la función son los dos modelos ajustados a comparar.
anova(res.cox8, res.cox7)## Analysis of Deviance Table
## Cox model: response is surv.obj
## Model 1: ~ +I(menostat == "Post") + I(NodeStatus == "NodePositive") + I(NodalStatus == "G2") + I(GrupoEdadCluster == "(50,66]") + I(Luminal == "TN")
## Model 2: ~ menostat + NodeStatus + NodalStatus + Luminal + GrupoEdadCluster
## loglik Chisq Df P(>|Chi|)
## 1 -4223.7
## 2 -4214.8 17.622 6 0.00725 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(res.cox7, res.cox6)## Analysis of Deviance Table
## Cox model: response is surv.obj
## Model 1: ~ menostat + NodeStatus + NodalStatus + Luminal + GrupoEdadCluster
## Model 2: ~ menostat + NodeStatus + Luminal + GrupoEdadCluster
## loglik Chisq Df P(>|Chi|)
## 1 -4214.8
## 2 -4218.7 7.6729 3 0.05328 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1En la tabla se muestra para los dos modelos la log-verosimilitud, el estadísticos Chisq, los grados de libertad Df y el p-valor P(>|Chi|) de la prueba.
3.5. Análisis exploratorio
3.5.1. Distrbución covariables con más significancia
Distribución de los eventos en el tiempo:
ggsurvevents(fit=survfit(surv.obj ~ 1, BC), normalized=T)
Las pacientes en las que tiene ocurrencia el evento son la inmensa mayoria, distribuidos estos uniformemente en el tiempo, con una diferencia de altura no demasiado grande para unos 25 meses.
Además, teniendo ya una previa idea de la importancia afecta cada covariable, haremos un análisis y lo veremos gráficamente.
Lo haremos de forma no paramétrica con Kaplain-Meier, ya que los demás no marcan diferencia respecto a este.
Empezaremos viendo más concretamente que antes la distribución de los datos marcados como más significativos (En %):
with(BC, prop.table(table(menostat, cens)))## cens
## menostat 0 1
## Post 0.08200734 0.67809058
## Pre 0.01468788 0.22521420with(BC, prop.table(table(NodalStatus, cens)))## cens
## NodalStatus 0 1
## G0 0.02692778 0.40881273
## G1 0.01468788 0.12239902
## G2 0.01223990 0.12239902
## G3 0.03182375 0.21052632
## G4 0.01101591 0.03916769with(BC, prop.table(table(Luminal, cens)))## cens
## Luminal 0 1
## HER2-enriched 0.006119951 0.035495716
## LuminalA 0.055079559 0.603427173
## LuminalB 0.014687882 0.119951040
## TN 0.020807834 0.144430845with(BC, table(Luminal, cens))## cens
## Luminal 0 1
## HER2-enriched 5 29
## LuminalA 45 493
## LuminalB 12 98
## TN 17 118with(BC, prop.table(table(NodeStatus, cens)))## cens
## NodeStatus 0 1
## Node negative 0.02692778 0.40881273
## Node positive 0.06976744 0.49449204with(BC, table(NodeStatus, cens))## cens
## NodeStatus 0 1
## Node negative 22 334
## Node positive 57 404En principio podemos ver que de las muertes acontecidas, las pacientes post menopausicas ocupan la gran mayoría, a su vez, en este mismo grupo los eventos no observados son muy pocos (un 0,008%).
También comprobamos lo que hemos visto en el modelo de Cox, además de Grado 0, las pacientes con G3 son aquellas en las que más muertes se han observado, y G4 en las que menos. De hecho si volvemos a los resultados de Cox, este nivel lo encontrabamos con NA.
La distribución de pacientes para estado nodel es más o menos homogenea, viendo unas pocas más de muertes en estado nodal positivo, aunque también son más las pacientes en esta categoría.
Solo con estas tablas no podemos decir demasiado, veámoslo con otro tipo de análisis.
Luminal tiene casi todas sus pacientes concentradas en LuminalA. El mayor número de muertes está en LuminalA y Triple negativo.
3.5.2. Análisis exploratorio no-paramétrico
survfit(surv.obj ~ menostat, BC, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia por menopausia", conf.int = T, legend.title = "Menopausia", pval=TRUE, ggtheme = theme_bw(), risk.table=T, risk.table.col="strata", risk.table.height = 0.25)
survfit(surv.obj ~ NodeStatus, BC, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia por estado nodal", conf.int = T, legend.title = "Estado Nodal", pval=TRUE, ggtheme = theme_bw(), palette = c("#E7B800", "#2E9FDF"), risk.table=T, risk.table.col="strata", risk.table.height = 0.25)
survfit(surv.obj ~ Luminal, BC, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia por grado", conf.int = T, legend.title = "Grado", pval=TRUE, ggtheme = theme_bw(), risk.table=T, risk.table.col="strata", risk.table.height = 0.32)
survfit(surv.obj ~ NodalStatus, BC, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia por grupo de nodos", conf.int = T, legend.title = "Estado de los nodos", pval=TRUE, ggtheme = theme_bw(), palette = "Dark2", risk.table=T, risk.table.col="strata", risk.table.height = 0.35)
Otra vez más vemos que el p-value de las diferencias entre menopausia marca la mayor diferencia entre esos grupos con notable diferencia al resto de variables, teniendo las pre mejor supervivencia.
En estado nodal la supervivencia para Nodo positivo es algo menor, pero no es tanta la diferencia, además, sus intervalos de confianza se solapan.
Con un pvalue también muy bajo pero no tanto como menopausia, encontramos el grupo del estado nodal. Aun con una buena diferencia según el pvalue, los intervalos de confianza son grandes y no podemos apreciar mucho, tal vez decir que al principio el grupo G4 y G3 tienen una supervivencia menor (como era de esperar) y los grupos G0, G1 y G2 mayor.
Para todas coincide que a partir del tiempo 100 quedan muy pocas pacientes vivas.
Con el menos fiable de todos los pvalues tenemos Luminal. Como habiamos visto. HER2-Herinched al tener tan pocos pacientes en ese grupo sus intervalos de confianza son muy amplios para poder alcanzar el 95% que le hemos marcado por defecto.
Para seguir con el estudio eliminaremos este grupo que causa problemas y dificultades.
BC<-BC[as.factor(BC$Luminal)!="HER2-enriched",]
BC$Luminal<-factor(BC$Luminal, labels=c("LuminalA", "LuminalB", "TN" ))
surv.obj<-Surv(BC$time, BC$cens)
survfit(surv.obj ~ Luminal, BC, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia por grado", conf.int = T, legend.title = "Grado", pval=TRUE, ggtheme = theme_bw(), risk.table=T, risk.table.col="strata", risk.table.height = 0.32)
El \(pvalue\) mejora, aunque no hemos conseguido que los intervalos dejen de estar solapados. LuminalA tiene la mayor parte de los valores, pero B y TN al tener menos sus intervalos son mayores. Por lo visto graficamente parece que la peor supervivencia es para LuminalB, y a contrario de lo que habiamos previsto, la mejor para TN. No son unas deduciones muy fiables ya que como vemos en la tabla de riesgo, que queden mmuchas menos pacientes vivas e debido a que en un principio habia menos que en LuminalA.
Las colas largas independientes que vemos en todos los gráficos, como bien se puede ver en las tabals de riesgo son debidas a que solo llega un paciente vivo a los últimos intervalos medidos. En el caso de grupo nodal solo sobrevive G2, para LuminalA, Nodalmente positivo y Post menopausica.
En cuanto a los grupos de edad:
survfit(surv.obj ~ GrupoEdad2, BC, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia por Grupos de edad", conf.int = T, legend.title = "Grupo edad", pval=TRUE, ggtheme = theme_bw())
survfit(surv.obj ~ GrupoEdadCluster, BC, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia por Grupos de edad Cluster", conf.int = T, legend.title = "Grupo edad Clustererizados", pval=TRUE, ggtheme = theme_bw())
Así como tgrade dejaremos de explorarlo que que demuestra no ser una aportación importante en el estudio, vemos que los grupos de edad si marcan mucha diferencia, siendo mejor y más fiables los clusterizados ya que sus intervalos son menores los agrupados que los hechos arbitrariamente. Los usaremos para seguir viendo diferencias.
Si aplicamos el modelo de cox para ver la significancia de todas las interacciones posibles:
summary(coxph(surv.obj~menostat*NodalStatus*Luminal, data=BC)) ## Call:
## coxph(formula = surv.obj ~ menostat * NodalStatus * Luminal,
## data = BC)
##
## n= 783, number of events= 709
##
## coef exp(coef) se(coef)
## menostatPre -0.27061 0.76291 0.17232
## NodalStatusG1 0.05886 1.06063 0.16337
## NodalStatusG2 0.09613 1.10090 0.15637
## NodalStatusG3 0.60381 1.82908 0.13382
## NodalStatusG4 0.21456 1.23931 0.34247
## LuminalLuminalB 0.44499 1.56047 0.20105
## LuminalTN 0.12469 1.13279 0.15610
## menostatPre:NodalStatusG1 0.14013 1.15043 0.31010
## menostatPre:NodalStatusG2 0.08851 1.09254 0.31403
## menostatPre:NodalStatusG3 -0.57836 0.56082 0.28163
## menostatPre:NodalStatusG4 0.59283 1.80911 0.46695
## menostatPre:LuminalLuminalB -0.38681 0.67922 0.48097
## menostatPre:LuminalTN -0.02697 0.97339 0.30969
## NodalStatusG1:LuminalLuminalB -0.45390 0.63515 0.41551
## NodalStatusG2:LuminalLuminalB -0.38070 0.68338 0.41260
## NodalStatusG3:LuminalLuminalB -0.86291 0.42193 0.30981
## NodalStatusG4:LuminalLuminalB -0.15493 0.85647 0.59333
## NodalStatusG1:LuminalTN -0.30393 0.73791 0.46074
## NodalStatusG2:LuminalTN 0.16604 1.18062 0.41010
## NodalStatusG3:LuminalTN -0.18090 0.83452 0.30305
## NodalStatusG4:LuminalTN -0.68125 0.50598 0.79733
## menostatPre:NodalStatusG1:LuminalLuminalB 0.79531 2.21512 0.74307
## menostatPre:NodalStatusG2:LuminalLuminalB -0.07554 0.92724 0.95515
## menostatPre:NodalStatusG3:LuminalLuminalB 1.28083 3.59962 0.67029
## menostatPre:NodalStatusG4:LuminalLuminalB NA NA 0.00000
## menostatPre:NodalStatusG1:LuminalTN -0.47980 0.61891 0.76252
## menostatPre:NodalStatusG2:LuminalTN -1.01026 0.36412 0.88979
## menostatPre:NodalStatusG3:LuminalTN 0.35075 1.42013 0.83802
## menostatPre:NodalStatusG4:LuminalTN 3.09207 22.02271 1.34340
## z Pr(>|z|)
## menostatPre -1.570 0.11633
## NodalStatusG1 0.360 0.71861
## NodalStatusG2 0.615 0.53870
## NodalStatusG3 4.512 6.42e-06 ***
## NodalStatusG4 0.626 0.53099
## LuminalLuminalB 2.213 0.02687 *
## LuminalTN 0.799 0.42443
## menostatPre:NodalStatusG1 0.452 0.65135
## menostatPre:NodalStatusG2 0.282 0.77807
## menostatPre:NodalStatusG3 -2.054 0.04001 *
## menostatPre:NodalStatusG4 1.270 0.20423
## menostatPre:LuminalLuminalB -0.804 0.42127
## menostatPre:LuminalTN -0.087 0.93061
## NodalStatusG1:LuminalLuminalB -1.092 0.27467
## NodalStatusG2:LuminalLuminalB -0.923 0.35617
## NodalStatusG3:LuminalLuminalB -2.785 0.00535 **
## NodalStatusG4:LuminalLuminalB -0.261 0.79400
## NodalStatusG1:LuminalTN -0.660 0.50948
## NodalStatusG2:LuminalTN 0.405 0.68557
## NodalStatusG3:LuminalTN -0.597 0.55056
## NodalStatusG4:LuminalTN -0.854 0.39288
## menostatPre:NodalStatusG1:LuminalLuminalB 1.070 0.28449
## menostatPre:NodalStatusG2:LuminalLuminalB -0.079 0.93696
## menostatPre:NodalStatusG3:LuminalLuminalB 1.911 0.05602 .
## menostatPre:NodalStatusG4:LuminalLuminalB NA NA
## menostatPre:NodalStatusG1:LuminalTN -0.629 0.52920
## menostatPre:NodalStatusG2:LuminalTN -1.135 0.25621
## menostatPre:NodalStatusG3:LuminalTN 0.419 0.67555
## menostatPre:NodalStatusG4:LuminalTN 2.302 0.02135 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95
## menostatPre 0.7629 1.31076 0.54424
## NodalStatusG1 1.0606 0.94283 0.77003
## NodalStatusG2 1.1009 0.90834 0.81031
## NodalStatusG3 1.8291 0.54672 1.40709
## NodalStatusG4 1.2393 0.80690 0.63339
## LuminalLuminalB 1.5605 0.64083 1.05226
## LuminalTN 1.1328 0.88277 0.83421
## menostatPre:NodalStatusG1 1.1504 0.86924 0.62647
## menostatPre:NodalStatusG2 1.0925 0.91530 0.59038
## menostatPre:NodalStatusG3 0.5608 1.78311 0.32292
## menostatPre:NodalStatusG4 1.8091 0.55276 0.72443
## menostatPre:LuminalLuminalB 0.6792 1.47227 0.26461
## menostatPre:LuminalTN 0.9734 1.02733 0.53050
## NodalStatusG1:LuminalLuminalB 0.6351 1.57444 0.28131
## NodalStatusG2:LuminalLuminalB 0.6834 1.46331 0.30441
## NodalStatusG3:LuminalLuminalB 0.4219 2.37004 0.22990
## NodalStatusG4:LuminalLuminalB 0.8565 1.16758 0.26771
## NodalStatusG1:LuminalTN 0.7379 1.35518 0.29910
## NodalStatusG2:LuminalTN 1.1806 0.84701 0.52848
## NodalStatusG3:LuminalTN 0.8345 1.19829 0.46077
## NodalStatusG4:LuminalTN 0.5060 1.97635 0.10603
## menostatPre:NodalStatusG1:LuminalLuminalB 2.2151 0.45144 0.51629
## menostatPre:NodalStatusG2:LuminalLuminalB 0.9272 1.07847 0.14262
## menostatPre:NodalStatusG3:LuminalLuminalB 3.5996 0.27781 0.96763
## menostatPre:NodalStatusG4:LuminalLuminalB NA NA NA
## menostatPre:NodalStatusG1:LuminalTN 0.6189 1.61576 0.13886
## menostatPre:NodalStatusG2:LuminalTN 0.3641 2.74631 0.06366
## menostatPre:NodalStatusG3:LuminalTN 1.4201 0.70416 0.27479
## menostatPre:NodalStatusG4:LuminalTN 22.0227 0.04541 1.58257
## upper .95
## menostatPre 1.0694
## NodalStatusG1 1.4609
## NodalStatusG2 1.4957
## NodalStatusG3 2.3776
## NodalStatusG4 2.4249
## LuminalLuminalB 2.3141
## LuminalTN 1.5382
## menostatPre:NodalStatusG1 2.1126
## menostatPre:NodalStatusG2 2.0218
## menostatPre:NodalStatusG3 0.9740
## menostatPre:NodalStatusG4 4.5179
## menostatPre:LuminalLuminalB 1.7435
## menostatPre:LuminalTN 1.7861
## NodalStatusG1:LuminalLuminalB 1.4340
## NodalStatusG2:LuminalLuminalB 1.5342
## NodalStatusG3:LuminalLuminalB 0.7744
## NodalStatusG4:LuminalLuminalB 2.7401
## NodalStatusG1:LuminalTN 1.8205
## NodalStatusG2:LuminalTN 2.6375
## NodalStatusG3:LuminalTN 1.5114
## NodalStatusG4:LuminalTN 2.4145
## menostatPre:NodalStatusG1:LuminalLuminalB 9.5039
## menostatPre:NodalStatusG2:LuminalLuminalB 6.0287
## menostatPre:NodalStatusG3:LuminalLuminalB 13.3908
## menostatPre:NodalStatusG4:LuminalLuminalB NA
## menostatPre:NodalStatusG1:LuminalTN 2.7586
## menostatPre:NodalStatusG2:LuminalTN 2.0828
## menostatPre:NodalStatusG3:LuminalTN 7.3393
## menostatPre:NodalStatusG4:LuminalTN 306.4633
##
## Concordance= 0.596 (se = 0.012 )
## Rsquare= 0.067 (max possible= 1 )
## Likelihood ratio test= 54.69 on 28 df, p=0.002
## Wald test = 59.92 on 28 df, p=4e-04
## Score (logrank) test = 69.53 on 28 df, p=2e-05Vemos que las combinaciones que mejor explcian el modelo son: menostat-tgrade, menostat-nodalStatus, tgrade-nodalStatus Haremos estas combinaciones teniendo como diferenciación principal el parámetro que ha demostrado ser el más significativo, menostat.
survfit(surv.obj ~ NodeStatus + menostat, BC, conf.type = "log-log") %>%
ggsurvplot(title = "Supervivencia estado nodal por menopausia", conf.int = T,
facet.by = "NodeStatus", legend.title = "Menostat", short.panel.labs = T)
survfit(surv.obj ~ NodalStatus + menostat, BC, conf.type = "log-log") %>%
ggsurvplot(title = "Supervivencia grupo nodal por menopausia", conf.int = T,
facet.by = "NodalStatus", legend.title = "Menostat", short.panel.labs = T)
survfit(surv.obj ~ Luminal+ menostat, BC, conf.type = "log-log") %>%
ggsurvplot(title = "Supervivencia Luminal por menopausia", conf.int=T,
facet.by = "Luminal", legend.title = "Menostat", short.panel.labs = T)
Y aplicando también la edad:
survfit(surv.obj~ NodeStatus+GrupoEdadCluster, BC, conf.type = "log-log") %>%
ggsurvplot(title = "Supervivencia de estado nodal por edad ", conf.int=T, palette = "Dark2",
facet.by = "GrupoEdadCluster")
survfit(surv.obj ~ NodalStatus+GrupoEdadCluster, BC, conf.type = "log-log") %>%
ggsurvplot(title = "Supervivencia entre Grupo nodal por edad", conf.int = F, palette = c("slategray1", "#2E9FDF","blue2", "#E7B800","orange3"), facet.by = "GrupoEdadCluster", legend.title = "Grupo nodal", short.panel.labs = T)
survfit(surv.obj ~ Luminal +GrupoEdadCluster, BC, conf.type = "log-log") %>%
ggsurvplot(title = "Supervivencia entre Luminal por edad", conf.int=F,
facet.by = "GrupoEdadCluster", legend.title = "Luminal", short.panel.labs = T)
3.5.2.1 Supuestos observados
Diferenciado por menopausia:
Cuando la paciente es nodalmente positiva, ser premenopausica marca una diferencia en supervivencia ligeramente mayor, pero a simple vista es bastante similar que post.
Para el grado de los nodos: Como ya hemos visto tantas veces y seguimos ocmprobando, G3 marca la mayor diferencia para los niveles de menopausia, siendo la supervivencia de las Pre mucho mayor. También hay una diferencia fiable en G0. Para G4 encontramos destacable que es el único donde la supervivencia es mejor para Post.
Para el Luminal: Como hemos visto en su correspondiente tabla de riesgo, Para LuminalB ,para el tiempo 100), solo quedan tres pacientes vivas que viendo la tabla de riesgo sabemos que es en el grupo Post, por eso ese corte tan extraño en Pre. En cuanto al resto siguen una distribución dentro de lo esperado, para los 3 la esperanza de vida es claramente menor en Post.
En general para esta covariable la esperanza de supervivencia para las Pacientes Premenstruales es mejor en la mayoría de los casos, exceptuando Grupo nodal G4. Es especialmente mejor en Grupos nodales G3 y G0 y Nodalmente positivo.
Siguiendo esta misma linea, esperamos que las pacientes del grupo de edad más jóvenes tengas mejores expectativas, coincidiendo con las premenopausicas.
Diferenciando por edad:
- Aquí encontramos la primera “contradicción”, así como veiamos que para las Pre la supervivencia era algo mejor con estado nodal positivo, viendo lor edad esto ocurre en el segundo grupo, donde todavia puede haber alguna menopausica, veamoslo:
with(BC, table(GrupoEdadCluster, menostat))## menostat
## GrupoEdadCluster Post Pre
## (25,50] 57 168
## (50,66] 328 20
## (66,90] 209 1Efectivamente el número de pre en el segundo grupo de edad es mucho menor, pero no por ello inexistente, por lo que no deja de tener sentido el supuesto que habiamos hecho.
El grupo que cuenta con más pacientes ((50,66]) tiene más probabilidades de sobrevivir siendo Node negative, y sus intervalos de confianza son menores, por lo que en principio nos fiaremos más de este supuesto.
Grupo nodal podemos ver el patrón de como la curva decae antes para los grupos G3 y G4 que para el resto, coincidendo en menor o mayor medida para los 3 grupos de edad.
Luminal no hay tanta diferencia como en las anteriores variables. Podemos destacar que Triple Negativo tiene mejor supervivencia en el primer Grupo de edad.
Todos los resultados obtenidos parecen coherentes y cercanos a lo que esperábamos.
3.6. Pruebas de hipótesis del análisis exploratorio
Tendremos dos hipótesis:
En este caso, dado que hemos visto que los intervalos de confianza son amplios ya que algunos niveles cuentan con pocos datos, estableceremos un nivel de confianza del \(99%\).
\(H_0\) o hipótesis nula, la cual no rechazaremos a no ser que se demuestre lo contrario (\(p-value \leq 0,01\) o probabilidad del 99% de que \(H_0\) sea verdadera): Independencia de las curvas de supervivencia.
\(H_1\) o hipótesis alternativa, la cual aceptaremos en caso de rechazar \(H_0\): Dependencia de curvas.
La función survdiff() nos permite realizar comparación por parejas de forma no paramétrica de funciones de supervivencia, (Es la forma interna que tiene survminer de calcualar el \(pvalue\) que veiamos representado en las gráficas). Podremos observar sus correspondientes \(pvalue\) y \(\chi^2\).
Como hemos visto anteriormente, podemos usar para el ajuste de LogRank
- rho=1: PetoPeto
- rho=0: logRank
Consiste en realizar un test estadístico para ver si los resultados se deben al azar.
LogRank nos da la similitud entre dos grupos. La suma de las dos últimas columnas nos darán el valor de \(\chi^2\).
Chi Square: \(\chi^2=\frac{(0_1-E_1)^2}{E_1}+\frac{(0_2-E_2)^2}{E_2}\), explica como son os resultados observados en comparación con los esperados. Es otra manera de aceptar o rechazar la hipótesis nula
Estimación de riesgo: \(OR=\frac{01/E1}{02/E2}\)
survdiff(surv.obj~ NodeStatus + strata(menostat), BC, rho=0)## Call:
## survdiff(formula = surv.obj ~ NodeStatus + strata(menostat),
## data = BC, rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## NodeStatus=Node negative 347 325 356 2.72 5.6
## NodeStatus=Node positive 436 384 353 2.75 5.6
##
## Chisq= 5.6 on 1 degrees of freedom, p= 0.02survdiff(surv.obj~ NodalStatus + strata(menostat), BC, rho=0)## Call:
## survdiff(formula = surv.obj ~ NodalStatus + strata(menostat),
## data = BC, rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## NodalStatus=G0 347 325 356.1 2.723 5.603
## NodalStatus=G1 108 96 102.8 0.451 0.534
## NodalStatus=G2 105 95 98.2 0.103 0.120
## NodalStatus=G3 185 163 129.5 8.678 10.894
## NodalStatus=G4 38 30 22.4 2.584 2.698
##
## Chisq= 15 on 4 degrees of freedom, p= 0.005survdiff(surv.obj~ Luminal + strata(menostat), BC, rho=0)## Call:
## survdiff(formula = surv.obj ~ Luminal + strata(menostat), data = BC,
## rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Luminal=LuminalA 538 493 501.9 0.1586 0.5480
## Luminal=LuminalB 110 98 86.2 1.6292 1.8800
## Luminal=TN 135 118 120.9 0.0707 0.0861
##
## Chisq= 1.9 on 2 degrees of freedom, p= 0.4survdiff(surv.obj~ NodeStatus+strata(GrupoEdadCluster), BC,rho=0)## Call:
## survdiff(formula = surv.obj ~ NodeStatus + strata(GrupoEdadCluster),
## data = BC, rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## NodeStatus=Node negative 347 325 355 2.58 5.32
## NodeStatus=Node positive 436 384 354 2.59 5.32
##
## Chisq= 5.3 on 1 degrees of freedom, p= 0.02survdiff(surv.obj~ NodalStatus+strata(GrupoEdadCluster), BC, rho=0)## Call:
## survdiff(formula = surv.obj ~ NodalStatus + strata(GrupoEdadCluster),
## data = BC, rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## NodalStatus=G0 347 325 355.3 2.5832 5.3188
## NodalStatus=G1 108 96 101.9 0.3455 0.4111
## NodalStatus=G2 105 95 97.5 0.0654 0.0768
## NodalStatus=G3 185 163 131.0 7.8125 9.8462
## NodalStatus=G4 38 30 23.2 1.9686 2.0617
##
## Chisq= 13.2 on 4 degrees of freedom, p= 0.01survdiff(surv.obj~ Luminal +strata(GrupoEdadCluster), BC, rho=0)## Call:
## survdiff(formula = surv.obj ~ Luminal + strata(GrupoEdadCluster),
## data = BC, rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Luminal=LuminalA 538 493 504 0.22650 0.7970
## Luminal=LuminalB 110 98 88 1.14202 1.3298
## Luminal=TN 135 118 117 0.00368 0.0045
##
## Chisq= 1.4 on 2 degrees of freedom, p= 0.5| Supuesto | Chisq | pvalue | Hipótesis |
|---|---|---|---|
| NodalStatus + menostat | 15 | 0.0005 | \(H_1\) |
| NodalStatus+GrupoEdadCluster | 13.2 | 0.01 | \(H_1\) |
| NodeStatus + menostat | 5.6 | 0.02 | \(H_0\) |
| NodeStatus+GrupoEdadCluster | 5.3 | 0.02 | \(H_0\) |
| Luminal + menostat | 1.9 | 0.4 | \(H_0\) |
| Luminal +GrupoEdadCluster | 1.4 | 0.5 | \(H_0\) |
Aceptamos \(H_0\) (independencia de las curvas) para un 99% de confianza para todas excepto para las estratificacionesmenopausia y edad para estado nodal, las cuales han rechazalo la hipótesis nula por tener un \(p-value<0,01\) y aceptamos su dependencia, entendiendo que no dan tanta información sobre el riesgo.
Parece que la estratificación por grupo de edad marca más diferencia que la menopausia, son grupos más específicos.
Nos quedamos con que la mejor explicación del riesgo viene dada por Luminal y Estado nodal (Positivo/Negativo), todos estratificados por los grupos de edad obtenidos de la clusterización.
survfit(surv.obj ~ Luminal +GrupoEdadCluster, BC, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia entre grado tumoral por edad", conf.int=F, facet.by = "GrupoEdadCluster", legend.title = "Grupo nodal", short.panel.labs = T, fun="cumhaz")
survfit(surv.obj~ NodeStatus+GrupoEdadCluster, BC, conf.type = "log-log") %>% ggsurvplot(title = "Supervivencia de estado nodal por edad ", conf.int=F, palette = "Dark2", facet.by = "GrupoEdadCluster", fun="cumhaz")
3.7. Análisis de estratificación de edad
3.7.1. Parametricamente vs No parametricamente
levels(BC$GrupoEdadCluster)## [1] "(25,50]" "(50,66]" "(66,90]"grupoedad1<-BC %>% filter(GrupoEdadCluster=="(25,50]")
grupoedad2<-BC %>% filter(GrupoEdadCluster=="(50,66]")
grupoedad3<-BC %>% filter(GrupoEdadCluster=="(50,66]")
grupos=list(grupoedad1,grupoedad2,grupoedad3)
params<- lapply(grupos, function(x) flexsurvreg(surv.obj ~1, data=x, dist = "exp") %>% summary(type="survival") %>% data.frame)
no.params<-lapply(grupos, function(x) survfit(surv.obj~ 1, data=x) %>% fortify)
gglist<-list()
gglist.risk<-list()
titles=c("(25,50]","(50,66]","(50,66]")
for(model in 1:length(grupos)) {
gglist[[model]]<-ggplot() + geom_line(aes(time, est, col = "Parametrico"), data=params[[model]]) + geom_step(aes(time, surv, col = "No Parametrico"), data=no.params[[model]]) + labs(x = "Tiempo", y = "Probabilidad de Supervivencia", col = "Ajustes", title = titles[model])
gglist.risk[[model]]<-ggplot() + geom_line(aes(time, -log(est), col = "Parametrico"), data=params[[model]]) + geom_step(aes(time, -log(surv), col = "No Parametrico"), data=no.params[[model]]) + labs(x = "Tiempo", y = "Riesgo acumulado", col = "Ajustes", title = titles[model])
}
library(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combinegrid.arrange(gglist[[1]], gglist[[2]],gglist[[3]], ncol=2, nrow=2)
grid.arrange(gglist.risk[[1]], gglist.risk[[2]],gglist.risk[[3]], ncol=2, nrow=2)
3.7.2. Estadísticamente
Los parámetros AIC, BIC y LogLik ya están explicados en el punto \(3.3.\)
params<- sapply(grupos, function(x) flexsurvreg(surv.obj ~1, data=x, dist = "exp"), USE.NAMES = T, simplify = F)
no.params<- sapply(grupos, function(x) survfit(surv.obj ~1, data=x), USE.NAMES = T, simplify = F)
means.param <- sapply(params, function(x) print(x, print.rmean=T))## Call:
## flexsurvreg(formula = surv.obj ~ 1, data = x, dist = "exp")
##
## Estimates:
## est L95% U95% se
## rate 0.024097 0.022387 0.025937 0.000905
##
## N = 783, Events: 709, Censored: 74
## Total time at risk: 29422.92
## Log-likelihood = -3350.503, df = 1
## AIC = 6703.005
##
## Call:
## flexsurvreg(formula = surv.obj ~ 1, data = x, dist = "exp")
##
## Estimates:
## est L95% U95% se
## rate 0.024097 0.022387 0.025937 0.000905
##
## N = 783, Events: 709, Censored: 74
## Total time at risk: 29422.92
## Log-likelihood = -3350.503, df = 1
## AIC = 6703.005
##
## Call:
## flexsurvreg(formula = surv.obj ~ 1, data = x, dist = "exp")
##
## Estimates:
## est L95% U95% se
## rate 0.024097 0.022387 0.025937 0.000905
##
## N = 783, Events: 709, Censored: 74
## Total time at risk: 29422.92
## Log-likelihood = -3350.503, df = 1
## AIC = 6703.005means.noparam<- sapply(no.params, function(x) print(x, print.rmean=T))## Call: survfit(formula = surv.obj ~ 1, data = x)
##
## n events *rmean *se(rmean) median 0.95LCL
## 783.0 709.0 40.7 1.3 27.6 24.5
## 0.95UCL
## 31.0
## * restricted mean with upper limit = 217
## Call: survfit(formula = surv.obj ~ 1, data = x)
##
## n events *rmean *se(rmean) median 0.95LCL
## 783.0 709.0 40.7 1.3 27.6 24.5
## 0.95UCL
## 31.0
## * restricted mean with upper limit = 217
## Call: survfit(formula = surv.obj ~ 1, data = x)
##
## n events *rmean *se(rmean) median 0.95LCL
## 783.0 709.0 40.7 1.3 27.6 24.5
## 0.95UCL
## 31.0
## * restricted mean with upper limit = 217La diferencia no es apreciable en ningún parámetro.
3.7.3. Cox
Proporcionalidad de los riesgos
cox.edad<-coxph(surv.obj ~ GrupoEdadCluster, BC)
check.ph<-cox.zph(cox.edad)
check.ph## rho chisq p
## GrupoEdadCluster(50,66] -0.0189 0.251 0.6162
## GrupoEdadCluster(66,90] -0.0747 3.862 0.0494
## GLOBAL NA 4.276 0.1179ggcoxzph(cox.zph(cox.edad))
Volvemos a comprobar que la hipótesis de riesgos proporcionales es válida para la estratificación de edad, y visualizamos la distrbución.
Schoenfeld residuals:
plot(check.ph, var = 1)
abline(h = coef(cox.edad)[1], col = "red", lwd = 2)
plot(check.ph, var = 1)
abline(h = coef(cox.edad)[2], col = "red", lwd = 2)
Vemos los 3 grupos claramente diferenciados, solo un de ellos encajando.
Análisis Residuales
ggcoxdiagnostics(cox.edad, type = "dfbeta")
Podemos ver que bastantes punyos tienen mucha influencia en la linea de regresión, y también tenemos algunos outlayers (más para el grupo de edad mayor)
ggcoxdiagnostics(res.cox, type = "deviance", ox.scale = "time")